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ABSTRACT 


The  desirable  performance  goal  of  an  ultrasonic  nondestructive 
evaluation,  NDE,  methodology  is  to  reliably  and  rapidly  obtain 
information  regarding  flaws  in  the  material  being  tested.  Decisions 
concerning  acceptance/rejection  of  material  for  further  usage  can  then 
be  made  on  the  basis  of  the  nature  and  severity  of  flaws  within  it. 
Flaw  size  estimates  are  currently  made  on  the  basis  of  an  idealized 
physical  model,  describing  the  flaw,  which  uses  the  computed  impulse 
response  as  an  input. 

The  research  reported  in  this  study  seeks  to  define  the  limits 
and  sensitivities  of  currently  available  deconvolution  algorithms.  In 
particular,  the  interrelationship  among  parameters  of  deconvolution 
procedures,  noise,  transducer  bandwidth  and  instrumentation  related 
errors  were  studied.  Simulation  experiments  dealing  with  the  impulse 
train  recovery  from  material  samples  are  illustrated  in  an  algorithmic 
manner  with  the  aid  of  a  laboratory  minicomputer  system. 
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INTRODUCTION 


Ultrasound  has  been  effectively  used  to  characterize  the  internal 
structure  of  objects  that  are  opaque  to  visible  light.  Historically 
the  information  contained  in  the  ultrasonic  signals ,  after  passage 
through  a  test  medium,  was  used  to  detect  the  presence  of  defects. 
Current  research  in  ultrasonic  testing  towards  a  nondestructive 
evaluation  capability  requires  extraction  of  information  concerning 
size  and  shape  of  flaws  from  experimental  data.  This  requires 
parametrization  of  the  defect  impulse  response  in  terms  of  physical 
scattering  mechanisms.  While  the  impulse  response  of  some 
geometrically  well  defined  flaws  (e.g.  a  spherical  flaw)  can  be 
theoretically  predicted,  a  theoretical  frame  work  for  arbitrary  shaped 
flaws  does  not  exist.  Recent  work  using  eigen-function  expansion 
(also  known  as  the  singularity  expansion  method)  for  canonical  objects 
such  as  spheres  or  spheriods  suggests  that  size  information  can  be 
extracted  from  the  flaw  impulse  response  in  a  pattern  recognition 
approach.  This  approach  is  akin  to  assuming  a  transfer  function  model 
and  extracting  exponential  terms  from  the  computed  impulse  response. 
In  both  of  these  methods,  described  in  the  literature  for  defect 
sizing  and  classification,  computation  of  impulse  response  becomes  a 
necessary  first  step. 

Modern  signal  processing  techniques  used  in  defining  flaw  impulse 
response  are  principally  based  on  comparative  analysis  of  a  reference 
signal  with  the  scattered  flaw  signal.  The  underlying  hypothesis  in 
these  techniques  is  that,  in  the  noise  free  case,  the  scattered  flaw 
signal  is  due  to  the  linear  convolution  of  the  ultrasonic  reference 
signal  with  the  flaw  impulse  response.  The  ultrasonic  reference 


signal  is  assumed  to  be  a  convolution  of  an  electrical  pulse  with  the 
instrumentation  impulse  response.  The  flaw  characterization  problem 
thus  reduces  to  determining  the  kernel  of  the  convolution  integral 
given  the  input  and  output  time  signals.  The  impulse  response 
recovery  (system  identification)  has  been  carried  out  both  in  the 
frequency  domain  (Wiener  filtering  approach)  and  recently  in  time 
domain  at  AFWAL/MLLP  (spline  function  approach).  One  of  the  authors. 
Dr.  P.  K.  Bhagat,  of  this  report  has  also  implemented  a  nonlinear 
processing  approach,  namely  homomorphic  processing,  at  the  Materials 
Laboratory  (AFWAL/MLLP)  for  impulse  recovery.  This  approach, 
utilizing  logarithmic  properties,  is  successful  in  identifying  both 
the  ultrasonic  pulse  used  as  well  as  the  sample  impulse  response  from 
simulated  backscattered  signal  alone.  Comparative  analysis  of  these 
three  techniques  with  synthetic  data  at  AFWAL/MLLP  reveals  that,  the 
recovered  impulse  responses  are  essentially  equivalent  in  defining  the 
time  separation  of  reflectors,  in  the  noise  free  case. 

One  theme  that  is  common  to  the  convolutional  model  is  that  the 
test  sample  is  assumed  to  be  a  lossless  medium,  which  is  not 
realistic.  Spectra  of  scattered  signals  from  actual  samples  indicate 
that  not  only  is  the  amplitude  of  the  received  signal  decreased,  but 
also  that  the  center  frequency  is  shifted.  The  shift  in  center 
frequency  cannot  be  explained  on  the  basis  of  a  linear  lumped 
parameter  lossless  model.  This  complicates  the  spectral  division 
processing  used  to  recover  the  flaw  impulse  response,  since  the 
dependence  of  the  interrogating  pulse  on  depth  in  the  medium  is  not 
accounted  for.  Equivalently,  in  the  time  domain  case,  changes  in  the 


It  has  been  suggested  that  the  test  sample  response  be  divided 
into  several  segments  and  the  impulse  response  recovered  on  the  basis 
that  the  interrogating  pulse  is  invariant  for  the  given  segment. 
While  this  model  is  certainly  more  realistic,  it  compounds  the  impulse 
response  recovery  problem  since  ultrasonic  pulses  appearing  in  the 
intermediate  segments  cannot  be  defined  experimentally. 

A  deconvolution  procedure  that  can  be  applied  to  recovery  of 
either  of  the  two  "short-time"  convolved  components  when  neither  is 
known  is  homomorphic  transformation.  The  assumption  required  under 
this  approach  is  that  the  complex  cepstra  for  impulse  response  and 
ultrasonic  pulse  do  not  overlap  in  the  cepstral  domain.  The 
ultrasonic  pulses  used  in  NDE  applications  have  a  smooth  spectrum  and 
tend  to  occupy  cepstral  space  around  the  origin  whereas  the  flaw 
impulse  response  function  appears  as  an  impulse  train.  The  spacing  of 
the  first  impulse  from  cepstral  origin,  known  as  first  arrival, 
determines  whether  the  cepstra  can  be  separated  or  not.  If  the  first 
arrival  and  complex  cepstrum  of  the  ultrasonic  pulse  do  not  overlap 
then  simple  gating  in  the  cepstral  space  can  be  used  to  recover  both 
the  impulse  response  and  the  ultrasonic  pulse  propagating  through  the 
medium. 

The  actual  signals  generated  as  a  result  of  measurements  in  NDE 
applications  are  never  noise  free.  A  more  realistic  model  would 
provide  for  convolutional  component  of  noise  due  to  coherent  clutter 
and  grain  scattering  plus  an  additive  component  due  to  electronic 
noise  and  random  experimental  errors.  Since  the  time/frequency/ 
cepstral  domain  techniques  are  based  on  a  linear  convolutional 
relationship  between  the  ultrasonic  pulse  and  sample  impulse  response 


straightforward  application  of  these  to  system  identification  will 
lead  to  errors  in  estimating  the  actual  impulse  response.  This 
problem  is  equivalently  known  as  ill-posedness  of  the  integral 
equation  in  the  time  domain  case  or  ill-conditioning  of  spectral 
division  in  the  frequency  domain.  In  the  cepstral  domain  the 
presence  of  convolutional  and  additive  noise  tends  to  blur  the 
boundary  between  the  pulse  and  impulse  response  cepstra  leading  to 
operator  dependent  recovery.  Thus,  in  all  the  deconvolution 
procedures  used,  the  impulse  response  recovery  tends  to  be  nonunique 
as  the  results  tend  to  depend  on  the  data  processing  requirements  of 
data  smoothing  and  filtering. 

While  there  has  been  a  major  thrust  in  ultrasonic  NDE  research 
towards  the  development  of  inversion  techniques,  very  little  published 
literature  exists  on  the  dependence  of  the  impulse  response  recovery 
on  the  signal  processing  procedures  used.  For  example,  the 
theoretical  impulse  response  of  a  weakly  scattering  ellipsoidal  object 
is  well  described  in  the  literature.  It  is  a  function  having  a 
constant  positive  amplitude  over  the  region  corresponding  to  the 
interior  of  the  object  and  two  negative  going  peaks  of  short  duration 
at  each  edge.  This  response  has  not  been  realized  even  for  test 
samples  containing  spherical  flaws  and  consequently  the  sizing 
algorithm,  "inverse  Born  approach",  has  led  to  erroneous  flaw  size 
estimates.  In  addition,  actual  flaws  may  have  rough  surfaces  which 
will  tend  to  broaden  the  peaks  in  the  recovered  impulse  response.  The 
deconvolution  procedure  used,  depending  on  the  choice  of  algorithm 
parameters,  will  introduce  additional  broadening  effects  which  will  be 
data-dependent .  Ever  present  noise  in  experimental  data  will  further 
degrade  the  impulse  response  recovery  process  depending  on  the 


particular  algorithm  used. 


It  is  apparent  from  the  foregoings  that  there  is  a  need  for 
development  of  a  practical  deconvolution  procedure  which  provides 
impulse  response  recovery  and  identification  based  on  a  realistic 
model  accounting  for  losses  in  the  medium  as  well  as  the  effects  of 
noise.  This  study  should  be  aimed  at  defining  the  limits  and 
character  of  convolutional  and  additive  noise  components  as  well  as 
segmentation  of  flaw  signal  for  optimum  impulse  response  recovery. 

In  an  attempt  to  define  the  range  of  applicability  and  practical 
implementation  of  the  current  deconvolution  procedures  for  NDE 
applications  a  simulation  study  was  carried  out  as  reported  here.  An 
integral  feature  of  this  study  was  the  close  collaboration  between  the 
principal  investigator  and  the  group  led  by  Dr.  T.  J.  Moran  at 
AFWAL/MLLP  in  the  development  and  implementation  of  signal  processing 
techniques  for  impulse  response  recovery  at  their  respective 
laboratories.  The  originally  proposed  three  year  study  was,  however, 
concluded  after  the  first  year  when  the  principal  investigator  (P.  K. 
Bhagat)  accepted  a  new  assignment  outside  the  University. 

The  overall  goal  of  the  research  study  detailed  herein  was  to 
establish  both  theoretically  and  in  a  practical  setting  the  strengths 
and  limitations  of  deconvolution  procedures  used  in  ultrasonic  NDE 
with  respect  to  noise  characteristics  and  instrumentation  employed. 
An  integral  feature  of  this  study  was  the  fact  that  the  procedures  are 
to  be  assessed  using  samples  containing  known  flaws  and  samples  where 
the  flaw  characteristics  are  unknown.  We  hope  that  the  approach 
outlined  here  will  lead  to  the  development  of  a  practical 
deconvolution  procedure  for  routine  field  inspections. 
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BACKGROUND 


The  earliest  approach  to  flaw  characterization  dates  back  to  1930 
where  the  amplitude  (intensity)  of  ultrasound  was  measured  after  its 
propagation  through  the  material  under  test.  Reduction  in  amplitude 
of  received  signal  compared  to  a  reference  was  interpreted  as  being 
caused  by  the  flaw.  Firestone  (1/2)  is  credited  as  being  the  first  to 
recognize  the  importance  of  the  pulse  echo  method  for  application  to 
nondestructive  evaluation,  NDE.  The  location  of  a  flaw  was  defined 
through  analysis  of  the  received  echo  pattern;  transit  time 
measurement  defines  flaw  location,  amplitude  of  the  echo  is  a  function 
of  flaw  size,  position  and  form  of  the  flaw  with  respect  to  the 
transducer  and  the  characteristics  of  the  instrumentation  used. 

In  general,  flaw  analysis  approaches  can  be  described  as  either 
parametric  or  nonparametric  methods.  By  parametric  methods,  we  mean 
those  techniques  which  involve  measurements  of  deterministric  physical 
parameters  resulting  from  the  material-sound  interaction  equations. 
Velocity  of  sound  propagation  and  attenuation  are  examples  of  these 
variables.  Nonparametric  methods  involve  essentially  statistical 
measurement  of  variables  which  generally  do  not  have  well  defined 
relationships  in  terms  of  physical  phenomena.  Feature  extraction, 
using  a  pattern  recognition  approach,  is  an  example  of  this 
methodology  which  has  been  used  in  the  author's  laboratory  for 
characterization  of  tissue  pathology  (3) . 

Quantitative  ultrasonic  NDE  techniques  are  based  primarily  on 
evaluation  of  impulse  response  (parametric  approach)  for  definition  of 
flaw  location,  shape  and  size.  As  has  been  described  earlier  the 
scattered  flaw  signal  Sf(t),  is  primarily  a  convolution  of  the 


interrogating  pulse,  sr(t),  with  the  flaw  impulse  response,  mf(t) 


sf(t)  =  sr(t)  *  mf(t)  (1) 

where  *  denotes  convolution 

Impulse  response  recovery  from  equation  1,  known  as  deconvolution,  has 

been  addressed  to  by  many  authors  and  appears  in  several  disciplines 
(for  example  see  references  4-13) .  In  the  frequency  domain  the  flaw 
transfer  function,  Mftjw),  may  be  written  as 

Mf(jw)  =  Sf(jw)  1 

Sr(jw)  (2) 

where  Sf(jw)  and  Sr(jw)  are  Fourier  transforms  of 

sf(t)  and  sr(t)  respectively  and  w  is  the  angular 

frequency. 

As  mentioned  earlier,  deconvolution  has  been  carried  out  in  the  time/ 
f requency/cepstral  domain  using  interactive  minicomputers.  In  the 
time  domain,  the  impulse  response  recovery  reduces  to  determining  the 
kernel  of  an  integral  equation  of  the  first  kind.  Phillips  (5)  has 
pointed  out  the  ill-posed  nature  of  this  problem.  Since  the  integral 
operator  does  not  have  a  bounded  inverse,  a  slight  change  in  measured 
data  will  cause  finite  changes  in  the  computed  kernel  values.  Thus, 
in  presence  of  noise,  deconvolution  using  equation  1  will  not  provide 
a  unique  solution.  Several  authors  (5,6)  have  developed  computer 
based  matrix  algorithms  which  perform  a  given  degree  of  smoothing  on 
the  data  to  minimize  the  effects  of  noise.  Strand  and  Westwater  (14) 
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have  developed  a  general  least  square  process  of  estimating  the 


solution  which  also  provides  a  measure  of  error  in  final  results. 


Hunt  (15)  has  developed  a  constrained  deconvolution  procedure  using 


discrete  Fourier  transformation,  which  is  equivalent  to  the  works  of 


Phillips  (5)  and  Twomey  (6)  .  Recently  Lee  (16)  has  proposed  a  two 


stage  solution  procedure  to  the  impulse  recovery  problem.  In  the 


first  step  the  reference  signal  is  fitted  in  the  least  square  sense 


using  spline  functions.  This  fitted  function  is  then  used  to  provide 


a  solution  to  the  deconvolution  problem.  Thus  the  user  can 


interactively  define  the  degree  of  smoothing  needed  for  his  data 


analysis  problem.  The  developed  algorithm,  due  to  the  choice  of 


spline  function,  is  quite  efficient  in  matrix  manipulations  and  has 


been  implemented  on  the  AFWAL/MLLP  computer  (35) . 


In  the  frequency  domain  the  problem  of  flaw  transfer  function 


recovery  reduces  to  the  design  of  an  inverse  filter.  Presence  of 


measurement  noise  ill-conditions  the  spectral  division  process  as 


shown  below: 


Mf(jw)  =  _ 1 _  [Sf(jw)  +  N2(jw)]  (3) 

Sr(jw) 


where  N2(jw)  is  the  noise  spectra. 


The  ratio  N2 ( jw)/Sr ( jw)  can  completely  dominate  the  Mf(jw)  computation 


in  frequency  bands  of  low  SNR.  In  order  to  account  for  noise  in  the 


scattered  signal  Wiener  filtering  and  constrained  deconvolution  have 


been  used  for  the  impulse  recovery  problem  (10,26) .  Wiener  filtering 


is  based  on  the  minimum  integral  mean  square  error  whereas  the 


constrained  deconvolution  uses  a  smoothness  constraint  function.  Both 


approaches  approximate  equation  2  by: 


w 


•s 


**Jvi 
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Mf(jw)  =  S  (jw)  Sf(jw) 

_ I _  (4) 

Sr  (jw)2  +  C2 

where  C  is  a  user  defined  parameter  and  Sr  (jw) 
is  complex  conjugate  of  Sr(jw) 

Since  neither  the  properties  of  noise  signals  nor  their  energy 
content  are  known  apriori,  practical  implementation  in  NDE  application 
chooses  C2  so  as  to  eliminate  spectral  division  in  low  SNR  regions. 
This  approach  has  been  used  extensively  in  quantitative  NDE  work  (27) . 
Furgason  et.  al  (26)  have  also  reported  on  the  deconvolution 
processing  for  flaw  signatures  and  provide  a  constrained 
deconvolution  filter  which  is  based  on  earlier  work  of  Hunt  (15)  . 
These  authors  suggest  that  the  deconvolution  process  can  be  made 
adaptive  if  the  target  response  can  be  separated  into  known  and 
unknown  components.  They  provide  impulse  response  computed  from 
synthetic  reflector  series  in  support  of  their  hypothesis. 

Goebbels  et.  al  (28)  have  considered  the  problem  of  grain 
scattering  present  in  actual  backscattered  signals.  They  formulate 
grain  scattering  as  a  superposition  of  convolution  outputs  from  each 
scatterer  excited  at  the  grain  boundaries  inside  the  sound  beam  and 
integrated  over  the  path  length.  Their  measurement  model  for 
scattered  signal,  x(t),  can  be  described  by: 


x  ( t )  =  sr(t)  *  [raf(t)  +  n1(t)]  (5) 

where  n^(t)  is  the  coherent  noise  component 

due  to  grains  in  the  material  sample. 


These  authors  (28)  minimize  the  grain  scattering  contributions 
through  spatial  averaging  of  the  backscattered  signal.  This  approach 
is  based  on  the  plausible  hypothesis  that  small  variations  in  the 
probe  position  cause  significant  changes  in  the  grain  scattered 
component  but  leave  the  component  due  to  flaw  essentially  unaffected. 
Obviously  this  hypothesis  implies  that  grain  size  is  much  smaller  than 
reflector  dimensions  and  ultrasonic  wavelength  used,  and  that  the 
grains  are  randomly  distributed  in  sample  space. 

In  the  cepstral  domain,  the  defining  equation  corresponding  to 
equation  2  is 

F”^  [log  Sf  (jw)  l-F'-^llog  Mf  ( jw)  ]+F-1[log  Sr  ( jw)  ]  (6) 

where  F-1[*]  is  the  inverse  Fourier  transform. 

Under  the  assumption  that  the  cepstra  of  the  flaw  impulse 
response  and  the  interrogating  pulse  occupy  disjoint  spaces  in  the 
cepstral  domain  algorithms  have  been  developed  for  deconvolution  in 
speech  processing,  image  processing,  ultrasonic  tissue 
characterization  and  seismic  analysis.  Recently  the  author  has 
implemented  a  deconvolution  procedure  based  on  equation  6  on  the 
AFWAL/MLLP  computer  (36).  Simulation  results  obtained  to  date 
indicate  that  cepstral  deconvolution  procedure  yields  comparable 
results  to  the  well  established  Wiener  filter  and  time  deconvolution 
procedure. 

While  limited  literature  exists  on  the  topic  of  homomorphic 
processing  applied  to  backscattered  data  analysis,  there  is  a  wealth 
of  papers  describing  its  application  in  other  areas.  Published 


applications  include  removal  of  image  blurs  (17,18,19),  speech 
processing  (11,20),  seisomology  (7,21),  and  ECG  signal  analysis  (22). 

There  have  been  several  studies  relating  to  the  effects  of 
additive  noise  on  cepstral  processing  in  the  literature.  Studies 
dealing  with  real  cepstrum  are  due  to  Cole  (17),  Cannon  (18),  and 
Stockham  (19) .  These  authors  were  mainly  interested  in  arriving  at  a 
power  spectrum  estimation  of  one  of  the  convolved  components  such  that 
a  suitable  gating  filter  could  be  implemented  for  deconvolution. 
Their  methodology  involved  averaging  the  real  cepstra  of  many  segments 
of  the  convolved  data  contaminated  with  noise.  They  did  not  consider 
affects  of  noise  on  a  single  segment  of  data.  Kemerait  and  Childers 
(23)  have  also  studied  the  degradation  caused  by  noise  on  the  real 
cepstrum.  They  showed  that  real  cepstral  peaks,  corresponding  to 
reflector  locations,  are  decreased  in  the  presence  of  noise.  The 
effect  of  additive  noise  on  the  complex  cepstrum  was  less  pronounced 
at  SNR  of  20  db .  The  work  of  Hassab  et  al  (24)  indicates  that 
successful  echo  detection  is  dependent  on  the  bandwidths  of  signal  and 
noise.  A  study  was  recently  carried  out  in  the  author's  laboratory  to 
investigate  the  effects  of  additive  noise  on  recovery  of  convolutional 
components.  An  ultrasonic  reference  pulse  was  convolved  with  an 
arbitrary  impulse  train  to  produce  a  synthetic  backscattered  signal. 
Scaled  Gaussian  noise  was  added  to  the  backscattered  signal  to  produce 
the  desired  SNR  in  data.  The  complex  cepstrum  was  computed  using 
exponentially  weighted  data  to  remove  the  undesirable  effects  of 
spectral  zeroes  in  unwrapping  the  phase.  The  results  obtained,  using 
symmetrical  gating  in  the  cepstral  domain,  indicate  excellent  recovery 
of  the  impulse  train  for  SNR  as  low  as  10  db  (25)  .  This  author  is 
unaware  of  any  published  work  which  takes  into  account  the  effect  of 
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coherent  noise  on  convolutional  component  recovery  in  cepstral 
processing. 


A  more  realistic  model  of  backscattered  data  from  actual  samples 
should  contain  contributions  from  both  coherent  and  incoherent  noise 
components.  Such  a  model  has  been  defined  by  Elsley  et.  al  (27)  who 
studied  low  frequency  characteristics  of  flaws  in  ceramics.  Their 
model,  defining  the  received  signal,  y(t) ,  in  the  time  domain,  is 
given  by: 

y  (t)  ■  s  j.  ( t)  *  [mf(t)  +  n^(t)  ]  +  n2(t)  (7) 

where  n^(t)  =  noise  component  defining  coherent 
clutter  and  grain  scattering  and, 
n2(t)  -  Electronic  random  noise. 
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Considering  each  noise  source  separately  with  the  scattering 
signal  sf(t)  these  authors  derive  a  constrained  deconvolution  filter. 
The  form  of  filter  (equation  4,  this  proposal)  is  essentially  the  same 
for  either  case,  differing  only  in  the  choice  of  C2.  They  conclude 
that  choice  of  C2  as  a  constant  based  on  results  obtained  by  their 
colleagues  is  satisfactory  for  deconvolution  purposes.  Elsley  and 
Addison  (29)  have  recently  studied  the  accuracy  of  one  dimensional 
Born  estimation  of  flaw  radius  in  the  presence  of  noise.  The 
synthetic  waveforms  used  in  this  study  were  the  calculated  scattering 
from  a  spherical  void  to  which  scaled  Gaussian  noise  was  added.  This 
noise  was  colored  to  accentuate  Rayleigh  or  grain  scattering 
[(frequency)4  dependence].  SNR  was  defined  as  the  square  root  of  the 
ratio  of  energy  in  the  flaw  signal  to  the  mean  energy  in  colored 
noise.  Impulse  response  was  computed  from  a  constrained  deconvolution 
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filter  and  integrated  in  the  frequency  domain  to  yield  the 
characteristic  function  corresponding  to  the  flaw.  Radius  of  the  flaw 
was  then  estimated  by  the  ratio  area/peak  value  of  the  characteristic 
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function.  For  a  flaw  diameter  of  800  ^m  a  plot  of  estimated  radius 


versus  SNR  is  given.  For  zero  dB  SNR  the  estimate  is  shown  to  be  off 


by  +  20%,  however,  one  must  note  that  this  case  implies  very  strongly 


scattering  data  from  coherent  noise  sources.  In  actual  practice 


preprocessing  will  reduce  the  extent  of  this  noise.  The 


characteristic  function  used  in  Born  approximation  has  also  been 


defined  as  an  impediographic  model  [Jones  (30)]  or  as  a  Raylograph 


[Beretsky  and  co-workers,  (31,  32,  33)].  These  authors  compute  the 


characteristic  function  as  an  integral  of  the  calculated  impulse 


response  and  also  suggest  an  iterative  procedure  for  optimization  of 


amplitude  and  time  location  of  the  recovered  impulse  response.  Their 


approach  involves  choosing  a  suitable  band  pass  filter  and  combining 


its  response  with  the  spectra  of  the  initial  impulse  response  recovery 


to  obtain  a  better  estimate.  While  this  analysis  is  shown  in  their 


paper  to  provide  excellent  results,  this  author  is  uncertain  of  the 


practical  utility  of  this  approach  for  quantitative  NDE. 


Bollig  and  Langenberg  (34)  have  approached  the  ultrasonic  defect 


classification  problem  in  terms  of  transfer  function  modelling.  These 


authors'  attempt  to  peel  off  exponentials  corresponding  to 


singularities  of  the  transfer  function.  These  authors  state  "The 


crucial  point  in  the  SEM-parametrization  of  experimentally  obtained 


time  records  of  scattered  ultrasonic  signals  is  still  the  lack  of  an 


efficient  algorithm  to  extract  the  singularities  out  of  the  transient 


data."  They  cite  Prony's  algorithm  which  has  been  used  to  extract  the 


singularities  from  the  transient  response  but  state  that  the  method  is 
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quite  sensitive  to  noise.  Thesis  authors  suggest  that  the 
singularities  may  be  used  in  a  feature  extraction  scheme  from  the 
system  impulse  response  as  a  pattern  recognition  procedure  for  flaw 
sizing.  They  give  results  identifying  creeping  wave  singularities  and 
scattering  cross  section  of  spherical  models  for  synthetic  data. 


MATERIAL  AND  METHODS 


The  experimental  set  up  for  data  acquisition  related  to  NDE  tests 
is  schematically  shown  in  Figure  1.  All  the  simulation  experiments 
were  conducted  in  a  pulse  echo  made  using  a  5  MHz  center  frequency 
spherically  focussed  transducer.  As  shown  in  Figure  1  a  PDP  11/23 
minicomputer  with  512  Kb  memory  and  479  Mb  disk  storage  is  the  major 
element  in  the  data  acquisition  procedure.  All  functions  which 
include  positioning  of  the  transducer  over  material  samples, 
energizing  the  transducer  and  capturing  the  digitized  return  echoes 
were  carried  out  with  the  aid  of  the  computer. 

An  MP  203  pulser  (Metrotek,  Inc)  and  an  UTA-3  transducer  analyzer 
(Aerotech  Labs)  were  used  to  excite  the  transducer  with  an  appropriate 
electrical  pulse,  gate  and  capture  the  reference  and  returning  echoes 
from  interfaces.  The  transducer,  sample  hold  and  the  reflector  were 
housed  in  a  plexiglass  tank  filled  with  distilled  water.  In  this 
study,  both  aluminum  block  and  plexiglass  reflectors  were  used. 

Three  dimensional  movement  of  the  ultrasonic  transducer  was 
facilitated  through  development  of  a  positioning  system.  This  system 
consists  of  three  independent  lead  screw  slide  assemblies  driven  by 
stepper  motors.  Number  of  pulses  input  to  the  motors  is  controlled  by 
user  written  software  commands.  Independent  shaft  encoders  were 
implemented  to  verify  and  accurately  position  the  transducers  in  a 
feedback  arrangement.  The  current  system  provides  for  a  normal 
incidence  of  the  transducer  on  the  sample  with  some  degree  of  probe 
rotation  (5  ym  accuracy) . 

Signal  digitization  was  carried  out  using  a  high  speed  A/D 
convertor  (Biomation  8100)  under  the  control  of  PDP  11/23.  This  A/D 


Fig fl.  .  Data  acquisition  and  computer  set  up 


convertor  features  a  maximum  sampling  rate  of  100  MHz  and  an  interval 
memory  of  2  K  words.  Under  software  control  the  digitized  data  was 
transferred  to  the  mass  storage  unit  of  PDP  11/23.  A  CRT  monitor 
(Tektronix  406)  was  provided  to  be  able  to  view  the  digitized  wave 
form  as  the  data  was  being  transferred. 

Data  analysis  was  performed  on  an  offline  basis  using  the  PDP 
11/23  computer.  Software  developed  for  the  three  deconvolution 
procedures  (listed  in  Appendix)  was  used  in  user  interactive  fashion 
to  generate  simulation  results.  The  programs  and  subroutines  used  in 
this  study  were  generally  adapted  from  literature  to  be  compatible 
with  PDP  11/23.  In  the  following  a  brief  summary  of  each  procedure 
usage  is  given. 

Time  Domain  Deconvolution  (Spline  fitting) 

Figure  2  shows  a  block  diagram  of  the  Time  Domain  Deconvolution 
Process.  As  shown  a  user  needs  to  specify  several  parameters  prior  to 
computation.  These  are;  the  knot  ratio  (KR) ,  the  order  of  solution 
spline  (OS) ,  and  the  number  of  basic  spline  function  (#BSP) .  KR,  once 
specified,  is  fixed  for  both  data  (the  reference,  s(ti),  and  the 
convolved,  y(t^) ,  data),  while  OS  and  #BSP  can  be  varied  independently 
by  user. 

For  the  data  used  in  this  experiment,  KR  values  of  3  or  4,  and  OS 
values  of  3  to  6  provide  good  fit.  Maximum  #BSP  is  50  due  to  the 
limited  memory  space  in  the  computer. 

The  computation  of  the  basic  spline  functions  are  carried  out  in 
the  subroutine  BSP.  The  aim  in  fitting  the  reference  data  is  to  solve 
the  linear  equation  Ac=b,  by  choosing  c  that  minimizes  the  least 
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am  of  the  Time  Domain  Deconvolution  process 


square  difference  between  the  original  and  the  smoothed  reference 
data.  Matrix  inversion  to  compute  c  =A“^b  is  done  in  the  subroutine 
BWS  (the  program  listing  is  provided  in  the  Appendix)  .  Once  the  Cj 
are  obtained,  these  values  are  then  used  to  smooth  the  y(t^)  data  and 
extract  the  impulse  response  approximate  solution,  h(t^)  from  it. 

Constrained  Deconvolution  (Wiener  filtering) 

The  Constrained  Deconvolution  Process  is  shown  in  Figure  3. 
Utilizing  the  DFT  routine,  all  the  analysis  is  done  in  the  frequency 
domain.  The  percent  cutoff  is  a  user  supplied  information,  and  is 
relative  to  the  reference  signal.  It  should  be  noted  that  in  contrast 
to  implementation  of  a  statistical  Wiener  filter  a  preselected  noise 
floor  is  used  for  deconvolution  in  this  study  (This  is  current 
practice  in  NDE) . 

By  examining  how  much  noise  contaminated  in  the  convolved  data,  a 
user  specifies  how  much  smoothing  should  be  done.  This  smoothing  is 
performed  by  setting  the  amplitude  spectrum  of  the  reference  data, 
which  has  the  value  less  than  or  equal  to  the  relative  percent  cutoff 
to  zero.  Notice  that  by  setting  the  amplitude  spectrum  to  zero,  the 
noise  which  usually  occupy  the  low  and  high  frequency  will  be  leveled 
off. 

Homomorphic  Deconvolution 

Figure  4  shows  The  Homomorphic  Deconvolution  Process 
schematically.  Some  important  points  to  be  considered  in  this  process 
are  given  below. 

1.  Append  zeros 

Appending  zeros  provides  two  advantages.  By  increasing 


Fig.  4  Block  diagram  of  the  Homomorphic  Deconvolution  process 


the  number  of  samples  aliasing  due  to  computation  logarithm  and 
phase  unwrapping  errors  caused  by  linear  phase  components  are 
reduced . 

Exponential  Weighting 

The  method  of  exponential  weighting  was  first  introduced  by 
Schafer,  discussed  in  [21],  to  convert  a  mixed  phase  sequence 
into  a  minimum  phase  sequence.  This  is  accomplished  through 
multiplication  of  input  data  sequence  by  a  real  number  an,  where 
0<a<l  and  n  is  the  data  sequence  number.  Choice  of  the 
weighting,  a,  is  data  dependent. 

Phase  Unwrapping 

Since  the  arctangent  routine  in  the  computer  only  evaluates 
the  phase  of  modulo  2  *  ,  discontinuities  occur  in  the  phase 
curve.  The  computation  of  the  logarithm  of  the  data  requires 
that  the  phase  must  be  analytic  in  some  annular  region  of  the  z- 
plane  [21].  Therefore,  phase  unwrapping  is  necessary  to  avoid 
these  discontinuities.  By  adding  the  appropriate  multiple  of  2  * 
to  the  principal  values  of  the  phase,  unwrapping  is  achieved. 

Filtering 

Linear  filtering  is  applied  in  the  complex  cepstrum  to 
decompose  the  convolved  sequences.  This  is  done  by  setting  the 
complex  cepstrum  due  to  the  contribution  of  one  or  the  other 
signal  to  zero.  Two  type  of  filters  have  been  used,  they  are 
low-pass  and  comb  or  notch  filter.  Low-pass  filtering  is  used 


when  the  first  arrival  is  easily  detected,  and  is  done  by  setting 
the  complex  cepstrum  to  zero  beyond  the  first  arrival.  The  comb 
filtering  is  used  when  the  complex  cepstrum  is  badly  corrupted  by 
noise  that  the  first  arrival  cannot  be  detected. 

The  Inverse  Process 

Inverse  processing  performed  to  transform  the  complex 
cepstrum  into  its  time  domain.  Reinsertion  of  the  linear  phase 
is  done  in  this  process  to  obtain  the  original  phase  from  its 
shifted  version  caused  by  the  linear  phase  removal.  Exponential 
unweighting  is  applied  to  the  sequences  to  restore  effect  of  the 
exponential  weighting  done  in  the  beginning  of  the  process. 

Noise  and  DFT  Routines 

In  order  to  make  the  data  more  realistic,  simulated  random 
noise  is  added  to  the  convolved  data.  The  random  noise  is 
generated  by  the  program  ADDNOI,  by  specifying  the  statistical 
level  of  the  noise  to  the  data. 

Simulation  Experiments 
1.  Noise  Free  Case 

Front  surface  echo  was  used  as  the  reference  signal.  This 
signal  was  convolved  with  arbitrary  impulse  trains  to  generate 
simulated  flaw  data.  The  polarity  and  relative  amplitudes  of  the 
impulses  were  defined  to  portray  backscatte red  signals  from 
inclusions  or  voids.  The  separation  between  impulses  was  varied 
according  to  the  ultrasonic  wavelength  used  in  an  effort  to 
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define  resolution  of  the  techniques  used. 

The  synthesized  signs-/  Sf(t),  was  deconvolved  using 
developed  programs  for  constrained  deconvolution,  cepstral 
processing  the  spline  function  deconvolution.  A  figure  of  merit, 
n,  defined  as: 

n  -  (T  [mf(ti)  -  mf(ti)]2dt 
Jo 

where  m^ft^)  is  the  recovered  impulse 

was  used  to  characterize  the  success/failure  of  each 
deconvolution  technique  used.  The  parameters  of  the 
deconvolution  procedure  used  were  varied  to  obtain  optimum 
results.  This  implies  varying  the  knot-spacing,  order  of  spline 
and  number  of  splines  in  case  of  time  domain  deconvolution; 
transducer  cut-off  frequency,  smoothing  operator  and  size  of  FFT 
in  case  of  constrained  deconvolution  (also  called  Wiener  filter 
in  ultrasonic  NDE)  and  exponential  weighting,  FFT  size  and  linear 
filter/gate  employed  in  cepstral  processing. 

In  addition,  the  impulse  response  recovery  problem  was 
studied  with  impulses  lying  within  an  ultrasonic  wavelength  of 
each  other  in  order  to  define  the  limitation  of  the  deconvolution 
procedures  in  identifying  closely  spaced  reflectors. 

2.  Coherent  Noise  Case 

Coherent  noise  was  generated  using  a  standard  Gaussian  noise 
algorithm  already  in  place  on  our  computer.  This  noise  in  some 
cases  will  be  colored  by  (frequency)^  weighting  to  simulate 
Rayleigh  scattering.  The  impulse  response  recovery  as  a  function 
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of  noise  level  was  studied.  Limited  experiments  have  already 
indicated  that  good  impulse  response  recovery  is  possible  for  SNR 
as  low  as  6  db.  However ,  these  experiments  were  carried  out  with 
wider  separation  between  impulses;  and  the  interactive  aspects  of 
both  noise  and  pulse  separation  were  not  fully  explored. 


RESULTS  AND  DISCUSSION 


Figure  5  shows  the  time  domain  plot  of  the  ultrasonic  reference 
signal  used  in  the  study.  This  signal  was  obtained  as  the  front 
surface  echo,  digitized  at  50  MHz,  from  an  aluminum  block  reflector. 
The  transducer  nominal  frequency  was  5  MHz.  As  can  be  seen  the  actual 
transducer  frequency  is  5.20  MHz  with  significant  high  frequency 
response  up  to  10  MHz.  Figure  6  shows  the  complex  cepstrum  of  data 
obtained  as  a  result  of  convolving  the  reference  signal  with  varying 
interface  (impulse)  separations.  At  an  interface  separation  of  25  * 
(wavelength) ,  the  first  arrival  (seen  as  a  small  sharp  peak)  can  be 
easily  discerned.  In  contrast  when  the  pulse  separation  was  reduced 
to  1.5  X,  the  first  arrival  is  submerged  in  the  cepstrum  of  the 
reference  signal.  Signal  separation  in  such  situations  is  difficult. 
Impulse  response  recovery  as  a  function  of  exponential  weight  is  shown 
in  Figure  7.  The  simulated  signal  represents  a  spherical  flaw  within 
a  material  sample.  Increasing  the  exponential  weight  tends  to 
smoothen  the  recovered  impulse  response.  The  presence  of  high 
frequency  data  in  the  results  is  due  mainly  to  inverse  logarithm 
(exponentiation)  function.  Also  note  the  displacement  of  time  origin 
of  data  due  to  linear  phase  (time  delay) .  Figure  8  shows  impulse 
response  recoveries  for  the  three  deconvolution  methods  KR=3  and  0S=4 
shows  good  recovery  but  there  are  a  lot  of  oscillations  in  the  data. 
In  the  absence  of  'apriori'  knowledge  of  interface  separations  it 
would  be  difficult  to  identify  the  structures  in  the  recovered  impulse 
response.  Constrained  deconvolution  (Wiener  filter)  using  Fourier 
transformed  reference  and  convolved  data  provides  excellent  recovery 
as  indicated.  Homomorphic  deconvolution  with  gate  width  of  20  samples 
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and  exponential  weight  of  0.9999  also  provides  good  recovery.  The 


presence  of  small  amplitude  high  frequency  oscillations  is  a  cause  for 


concern.  At  a  SNR  of  30  dB  the  impulse  response  recovery  is  degraded 


as  shown  in  Figure  9.  Note  the  smoothing  function  of  the  time  domain 


deconvolution.  There  is  no  perceptible  change  in  the  shape  of  the 


recovered  impulse  train.  In  contrast,  both  the  Wiener  and  cepstral 


processed  data  show  significant  noise  floor  in  the  recovered  impulse 


train.  Figure  10  and  11  show  the  impulse  recoveries  for  20  dB  and  10 


dB  SNR.  Degradation  of  impulse  response  recovery  is  apparent  in  all 


the  methods.  Interestingly,  at  10  dB  SNR  time  domain  deconvolution 


provides  sharpened  peaks  coupled  with  low  frequency  interference.  In 


contrast,  homomorphic  deconvolution  shows  large  amplitude  noise  in  the 


recovered  data.  Table  1  shows  the  normalized  RMS  error  in  recovering 


the  impulse  train.  Figure  12  shows  the  percentage  deviation  of  Born 


estimate  of  flaw  dimensions  based  on  recovered  impulse  responses. 


Degraded  performance  is  apparent  in  all  the  three  deconvolution 


procedures  used. 


WWW 


-vXi ZnAr.t 


c 

U! 


Hi 


2  O  (9 

IU  fi  & 
—  I  Hi 

3  h  o 

4+0 


4  o  + 


+  04 


o  4  + 


o  4  4- 

X 


(0 

DC! 


# 

w 

> 

0) 

u 


a 

w 

a 

u 

o 

CD 


<4-1 

o 

c 

o 


> 

<u 

o 

S'* 


3 

3 


o 

C9 


O 

CM 


o 


o 

M 

6h 


SUMMARY 


An  NDE  test  setup  with  capability  of  computer  based  ultrasonic 
flaw  signal  analysis  is  detailed  in  this  study.  Extensive  software 
development  along  with  fabrication  of  a  three  dimensional  positioning 
system  formed  a  major  portion  of  the  study.  Limited  simulation 
results  indicate  that  the  deconvolution  procedures  are  greatly 
dependent  on  the  SNR  of  data.  Care  must  be  exercised  in  using  the 
deconvolved  data  to  estimate  flaw  dimensions  due  to  the  presence  of 
undesirable  noise  in  the  recovered  impulse  response.  Other  approaches 
such  as  autoregressive  modeling  of  the  signal  might  provide  a  vehicle 
for  noise  suppression  and  smoothing. 


SIGNIFICANCE  OF  RESEARCH 


Careful  documentation  of  sensitivities  and  limitations  of 
deconvolution  procedures  for  impulse  response  recovery  will  constitute 
a  major  contribution  to  the  field  of  quantitative  ultrasonic  NDE. 
Since  all  the  physical  models  developed  to  data  for  flaw  sizing 
require  computation  of  sample  transfer  function  as  a  necessary  first 
step,  the  research  outlined  here  will  provide  both  a  theoretical  basis 
as  well  a  practical  limitations  on  impulse  recovery  process.  If  we 
are  successful  in  developing  a  practical  deconvolution  procedure  for 
real  time  application  under  actual  field  conditions,  significant 
improvement  in  flaw  sizing  capability  will  result.  The 
interdisciplinary  nature  of  the  project  and  collaboration  with  Air 
Force  personnel  should  serve  to  promote  NDE  research  directly-  oriented 
towards  Air  Force  needs  and  increase  University  participation  in  this 
kind  of  activity. 
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APPENDIX 


Software  Listings 

Complete  FORTRAN  listings  of  the  deconvolution  procedures  used  in 
this  study  are  provided.  The  programs  are  generally  in  subroutine 
form  and  therefore  can  be  adapted  to  different  computer  systems  with  a 
minimal  amount  of  modifications. 


SUBROUTINE  CCEPS < NX , X t ISNX , I SFX > I SSUC t CX r AUX > 

U  _ 

U  The  subroutines  used  in  this  Program  were  taken  from: 

w 

'Proslrams  for  Digital  Signal'  Processing' 

C  IEEE  Press?  197? 

C  345  East  47  Strettj  New  York?  NY  10017 

C  Sponsored  by  the  IEEE  Acoustics!  Speech?  and 

C  Signal  Processing  Society 

C  Lib.  of  Congress  Cat,  Card  *  79-89028 

C  IEEE  Book  S:  0-37942-128-2  (paperback  ver.) 

C  #  0-37942-127-4  (hardback) 

C  Also  published  by  John  Wiley  S  Sons?  Inc. 

C  Wiley  Order  *  0-471-05961-7  (paperback  ver.) 

C  t  0-471-05962-5  (hardback) 

DIMENSION  X(l)fCX(l)f AUX ( 1 ) 

COMMON  PI?TU0PI?THLINC?THLC0N?NFFT?NPTS?N?L?H?H1?DVTMN2 
LOGICAL  I3SUC 
C 

C  SUBROUTINE  TO  COMPUTE  THE  COMPLEX  CEPSTRUM 

C 

NPT3=NFFT/2 
N  =  1 2 
L=2**N 

H=FLOAT(L)*FLOAT(NFFT) 

H 1 =PI /H 
ISSUC=.TRUE, 

ISNX=1 

C  TRANSFORM  X(N)  AND  N*X(N) 

C 

DO  10  1  =  1 ?  NX 
CX(I)=X(I) 

AUX(I)=FL0AT(I-1)*X(I) 

10  CONTINUE 

C  APPEND  THE  NECESSARY  ZEROES  FOR  FFT  COMPUTATION 
C 

INITL=NX+1 
IEND  =  NFFT +2 
DO  20  I=INITL? IEND 
CX(I>=0.0 
AUX ( I ) =0 , 0 
20  CONTINUE 

C 

C  COMPUTE  FFT 

CALL  FAST ( CX  ?  NFFT  > 

CALL  FAST ( AUX ?  NFFT ) 

C  CHECK  IF  SIGN  REVERSAL  IS  REQUIRED 
C 

IF ( CX ( 1 ) .LT.0.0) ISNXa-1 
C 

c 

c 

C  COMPUTE  MAGNITUDE  OF  SPECTRUM  5  STORE  IN  ODD  INDEXED  VALUES  OF 
C  AUX 
C 

C  COMPUTE  PHASE  DERIVATIVE  OF  THE  SPECTRUM  J  STORE  IN  EVEN  INDEXED 
C  VALUES  OF  AUX 
C 


I0  =  -1 


1 


on  n  n  o  o-  no  o  ( ri  i 


’  'v-V 


30 


c 

c 


c 

c 

c 


AO 


0 

1 

0 


DVTMN2=0 . 0 
IENB=NPTS+1 
no  30  I=1>IENB 
:u=i0f2 
IE-10+1 

AMAGSQ=AM0DSQ(CX(I0> ,CX(IE) ) 

PDVT*PHADVT<CX(IQ) >CX< IE) ,AUX( 10) r AUX (  I E  )  »  AMAGSQ > 

AUX< IQ)=AMAGSG 

AUX ( IE ) =PQVT 

BVTMN2*DVTMN2+PDVT 

CONTINUE 

DVTMH2*(2.*DVTMN2-AUX(2)-PDVT) /(FLOAT! NPTS) ) 


PPBVT =AUX ( 2 ) 

PPHASE=0 ♦ 0 

PPV=PPVPHA(CX(1 ) >CX<2) i I SNX ) 

CX(1)=0.5*AL0G(AUX(1) ) 

CX ( 2 ) =0 . 0 
10  =  1 

DO  50  I  =  2  >  I  £  M  D 
10  =  104-2 
I E  = 1 0  +  1 
PDVT=AUX( IE; 

PPV*PPVPHA(CX(IO> » CX ( IE ) tISNX) 

.  PHASE=PHAUNW(X»NX> ISNX* I »PPHASEf PPDVT » PPV> PDVT r ISSUE) 

IF  PHASE  ESTIMATE  SUCCESSFUL  ? CONTINUE 

IF  < ISSUC ) GOTO  40 
IS3UC=. FALSE. 

RETURN 

PPDVT  =  PDUT 

PPHASE=PHASE 

CX( I0)=0.5*ALOG(AUX( IQ) ) 

CX ( IE ) =PHASE 
TYPE  **IE»CX(IE) 

TYPE  41>'033»*133»*101> I FIX (FLOAT ( I > /FLOAT ( IEND>#100 ) 
FORMAT  (  '  '  »3A1 r 13  f '  V.' ) 

CONTINUE 


REMOVE  LINEAR  PHASE  COMPONENT 


ISFX=(ABS(PHASE/PI)+.l) 

IF  (PHASE. LT. 0. 0) ISFX=-ISFX 
H= PHASE/ FLO AT (NPTS) 

IE=0 

DO  60  1=1 » IEND 
IE=IE+2 

CX(IE)=CX(IE)-H*FLQAT(I-1) 

CONTINUE 


COMPUTE  THE  COMPLEX  CEPSTRUM 

CALL  FSST(CX>NFFT) 

RETURN 

END 


SUBROUTINE  3PCVAL( NX> X » FREQ » XR » XI » YR » YI > 


rv;;  •  *v  v-  v*  wy  \*y  vv v 


D t MENS  ION  X(i) 

DOUBLE  PRECISION  UO > U1 » U2 » WO» W1 » W2 » A» B» C > D » A1 » A2 » 
1 3 AO  7 CAO  .  X J 


I N I  T 

C*iOaDBLE(  COS  (FREQ)  > 

3 A 0  ~ 0 BLE (SI M ( F R  E  Q ) ) 

i'il  =2 .  D  +  0*CA0 

U 1  =  0  »  D  0 

U2  =  U1 

U1--U1 

l-J  2  =  U 1 

DO  10  J - 1 »  N X 
'< J  =  D  B  L.  E  <  X  <  J }  > 

U0-XJ  +  A1.HU1-U2 

i-IO=  CHBLE(  FLOAT  (  J-l  )  )  )  *X J  +  A 1  *U  1 -W2 
U2  =  U1 

m=uo 

W  2  ■=  W 1 
U1  =  W0 
CONTINUE 


A~U1-U2*CA0 
B=U2*SA0 
C~W1-W2*CA0 
D  W  2  $  3  A  0 

A2=DBLE ( FREQtfFLQAT ( NX-1 )  ) 
U1=DC0S ( A2 ) 

U2“-DSIN< A2> 
XR=SNGL<U1*A-U2*B> 

XI  =SNGL(U2*A-!-Ul*&> 
YR=SNGL<U1*C-U2*D> 

Y [=SNGL(U2*C+U1*D> 

RETURN 

END 


FUNCT  ION  F’HAUNW  <  X  r  NX  f ISNX » I > F’PHASE  »  F’PDVT  r  PPV  »  F’DVT  » I  SCONS  ) 

ROUTINE  TO  DO  PHASE  UNWRAPPING 

DIMENSION  SDVT  ( 17 )  »  SF'F'V <  1 7 )  »  X<  1 ) 

INTEGER  SINBEX(17)  .F’lNDEX.SP 
LOGICAL  I  SCONS. FIRST 

COMMON  P I  .  TWOP I .  THLINC  >  THLCON  t  NFFT  ;  NF’TS  »N»L»H»Hi»  DUTMN2 


FIRS T=. TRUE, 

PINDEX=1 

SP*1 

3PPV(SP)=PPV 
3  DVT  <  S’  P  >  =  P  D  V  T 


3 


.f  INDEX  (SP>sL+l 


c::j  ru  10 

!•  INDEX  =S  I  NDEX-:  SP  ) 
PPHASE=PHASE 
:;,:'UV  f  =SDVT<SP) 

;3  r  ~  3  P  -  1 
S(J  TO  TO 


C 

• ;  a 

—  v 


C 

30 

C 


c 


c 

>10 


cc 

cc 


r* 

w 


CC 

c 

c 


c 

c 

c 


IF ■ (SINDEX <  SP ) -P INDEX ) «  GT . 1 ) GO  TO  30 
ISCQNS=. FALSE. 

PHAUNU-O, 

RETURN 


i<-  <  SINDEX  (  SP )  +P INDEX )  /2 

FREQ=TWOPI* (FLOAT ( 1-2 ) *FLOAT < L > +FLOAT (K-l ) >/H 
CALL  SPCVAL(NX»X»FREQ»  XR»XI r YR» YI > 


SPaSP+1 
SINDEX  <  SP ) ~K 

SPPV(SP)=PPVPHA(XR»XI» ISNX ) 

X  M  A  G  =  A  M  0  D  S  Q  (  X  R  ?  X I ) 

SOVT ( SP ) =PHADVT ( XR »  XI »  YR »  YI » XMAG ) 

DELTA=H1 AFLOAT ( SINDEX ( SP ) -PI NDEX ) 
PHAINC=DELTA*(PPDVT+SDUT(SP)  ) 


IF( ABS(PHAINC-DELTA*DVTMN2) . GT , THLINC ) GOTO  20 
PHASE=PPHASE+PHAINC 

CALL  PHCHCK(PHASE»SPPV(SP) f ISCONS) 

IF  ( .NOT. I3C0NSJG0T0  20 


IF' ABS(PHASE-PPHASE) .GT.PI )G0  TO  20 


IF { SP ♦ NE . 1 ) GOTO  10 

PHAUNU=PHASE 

RETURN 

END 


FUNCTION  PPVPHA<XR»XI ? ISNX) 

IF ( XR ♦ EQ ♦ 0 . 0  .AND.  XI  ,EQ.  0,0)  PPVPHA=0.0 
IF ( XR . EQ . 0 . 0  .AND.  XI  .EQ.  0.0)  RETURN 

IF( I3NX.EQ. 1 )  PPVPHA=(ATAN2(XIfXR) ) 

IF ( ISNX . EQ . (-1) )  PPVPHA=(ATAN2(-(XI ) ,-<XR) ) 

RETURN 

END 

4 


o  o  o 


u 


FUNCTION  PHADVT'XRfXIfYRfYIfXMAG) 

IF ( XMAG  . EG ♦  0,0)  PHADVT=0.0 
IF ( XMAS  .EQ,  0.0)  RETURN 

PH ADVT=-SNGL( i DELE \ XR ) $DBLE C YR ) +DBLE < XI ) #DBLE  <  YI ) > /DBL£<  XMAG 

RETURN 

END 


) 


FUNCTION  AMODSQ  (ZRfZI) 

AMODSQ=SNGL<DBLE(ZR>*DBLE<ZR)+DBLE<ZI )*DBLE(ZI )  ) 

RETURN 

END 


SUBROUTINE  PHCHCK ( PH i PV » ISCONS ) 

C 

C  THIS  ROUTINE 

C  CHECKS  THE  PHASE  DIFFERENCE  BETWEEN  THE  PREVIOUS  DATA  POINTS  AND 

C  IF  THE  PHASE  DIFF  DOES  NOT  EXCEED  THE  USER  DEFINED  THRESHOLD 

C  THEN  THE  PHASE  VALUE  ISN'T  CHANGED. 

C 

COMMON  PIfTUOPI » THLINC f THLCON » NFFT t NPTS f N f L f H f HI fDVTMN2 
LOGICAL  ISCONS 

C 

AO=(PH-PV)/TUOPI 

A1=FL0AT ( IFIX(AO) )*TWOPI+PV 

A2=A1+SIGN<TU0PI >A0) 

A3=ABS< Al-PH) 

A4=ABS ( A2-PH ) 

C  CHECK  CONSISTENCY  . 

ISCONS*. FALSE. 

IF < A3. GT.THLCON. AND. A4.GT.THLC0N) RETURN 
T3CCNS= . TRUE ♦ 

r* 

P  H = A 1 

IF  < A3.GT,A4)PH  =  A2 

RETURN 

END 


SUBROUTINE  ICEPS ( CX f ISNX f ISFX > 

DIMENSION  CX ( 2 ) 

COMMON  PI f TWOPI f THLINC fTHLC0NfNFFTfNPTSfNfLfHfH1fDVTMN2 

CCCC 

CC 

c 

SNX=FLOAT (IS NX) 

SFX=H/2 . 

CX ( NFFT+1 ) =0 . 

CX ( NFFT  +  2 ) =0 , 

CALL  FAST  <  CX  »  NFFT ) 

CX<1)=SNX*EXP(CX(1> ) 

C 

C  PERFORM  EXPONENTIATION  OF  TRANSFORMED  DATA 


5 


J  u  1  j  ;  -  j  j  i  J  r  r  i  t  i  >  1' 
i  ~K  r  1 

P  D  L  Y  -  S  F  X F  L  u  A  T  (  K  - 1  ) 
r-5NX  *EXP- CX(!U  > 

C  A  <  K )  -  T t C 0 S  <  C X  •:  K  l  •  *  P  H  D  L  Y  ) 
c ;;  *  K 1  ‘  -T  x  S I H  <  C  X  <  K 1  >  +  P  H  D  I-  Y  > 

CONTINUE 

MLW  PERFORM  INVERSE  FOURIER  TRANSFORM 

CX  < NFF7  +  21 =0  , 

CALL  FSST ( CX »  NFFT ) 

RETURN 

END 


PROGRAM  TESTGT 

Tins  program  is  used  to  provide  convolved  date  for  use  in 
honomo rphic_ p rocess i ng , the  user  provides  an  ultrasonic  pulse  to 
be  used  as  reference  and  a  synthetic  reflector  series. the  output 
of  this  program  is  a  synthetic  flaw  signal ♦ synthetic  reflector 
series  is  defined  at  integer  sAmples  of  time. 

COMMON  PI » TUOPI » THLINC *  THLCON » NFFT r MPTS  »  N  f  L » H »  HI » DVTMN2 
DIMENSION  CX( 1026) »AUX( 1024) 

LOGICAL  1  FILNAM  (  15  )  »  TRUE  >  FALSE 
DATA  TRUE/. TRUE./  FALSE/ . FALSE . / 

CALL  JT I TLE (  '  TESTGT  'tot »2.0) 

GET  PULSE  DATA  FILE  NAME 

jvpE  25 

F0RMAT('  Enter  ultrasonic  reference  data  filename J  '?£) 

CALL  GETDFN  <  FILNAM ) 

IF  (0PNFIL( 3 » FILNAM » 'R' ) .NE.TRUE)  GOTO  20 
CALL  GETDAT(3fNPf TFCTRf CXfTRUE) 

READ  THE  REFLECTOR  SERIES  DATA 

TYPE  30 

F0RMAT< '  Enter  synthetic  reflector  series  data  filename;  '»$) 
CALL  GETDFN ( FILNAM ) 

IF  <  GPMFIL ( 2  » FILNAM  ? ' R ' ) • NE • TRUE )  GOTO  26 
CALL  GETDAT(2fNP1 f TFCTR1 » AUXf .TRUE. ) 

NP2  -=  NP 

IF  <  NPl . GT. NP2 )  NP2  =  NP1 
N P 2  =  NP2+NP2 
TYPE  35 

FORMAT*'  Enter  size  of  DFT» length  must  be  a  power  of 
NFFT  =  i PROMT ( NP2 ) 

Check  for  correct  NFFT 

IF { NFFT . GE . 4NP+NP1 ) >  GOTO  40 

TYPE  fr'  For  linear  convolution  dft  length  must  be 
1;  or  =  sum  of  dats  points  in  both  files!' 
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:  n  :  t  l  --  n  p  i  + 1 

DO  125  I=INITLfIEND 
A  U  X  (  I  >  ~  0 . 0 

COMPUTE  F FT 

CALL  FAST- CX»NFFT) 

CALL  FAST (AUXfNFFT) 

Commute  line 3 r  e o n v olution  .  Make  s u  p e  t h a t  t h s  d f  t  1  e n S  t  r 
i a  Greater  than  '  NPfNPl-1 '  and  a  lac  a  multiple  of  two. 

T Y  P  E  'Computing  1  i  n e  a  r  convolution...' 

DC  130  1  =  1? NFFT  +  1  >  2  !  This  loop  computes  U 

T1-CX( I ) #AUX  ( I )  -  CX  < 1  +  1 )*AUXv 1  +  1 )  !  product  o  f  two  D  F  T  '  s ! 

T2=CX(  1  +  1 >*AUX< I )tCX( I >*AUX< I  +  l >  !  CX  and  AUX 

CXv I ) =  T 1 
CX( I+l >=T2 
CONTINUE 

T  a k s  inverse  fourier  t  pans  fo r m  a  n d  write  t h e  c  o n v o 1 u  t i o n 
r  e suits  out. 

CALL  F3ST(CX»NFFT> 


U  p i t  e  t h e  data  for  plotting  purposes 


F G R M A T \ / »  Enter  output  c o n volved  ti m e  d a  t 
CALL  GETDFN(FILNAM) 

IF  (  GF'NFIL  (  4  f  FILNAM  >  '  W  '  >  .  NE  .  TRUE  )  GOTO  45 
C A LL  P U T D A T ( 4  r  NFFT r TFCTR  *  C X ) 


r  Kenanie. 


n  o  t2 


C 

c 

c 

w 

c 
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30 


jf\ 


PROGRAM  AD DM 01 

Tiiis  program  accepts  a  predefined  filename  '.done rail?  -jet j 
f  ro«i  a  reflector  wave  end  computes  its  statistics.  Based 
upon  these  statistics  and  j  SMR  value  input  be  the  user* 
noise  data  is  produced.  The  statistics  of  the  noise  data 
are  then  computed  including  the  actual  SNR  ratio.  The  noise 
data  is  Gausisn  distributed  using  central  limit  theorem. 

-.Gregg  Liming  3/*/35 


(REF  I- 1 1,  DAT) 

(REF***  .N;::;) 
where  x;;  =  SNR 
( IMP***. DAT) 

(  IMP***.Nxx> 
where  x x  -  SNR 

Document ! * document !  What  would  Dr  B.  say!  3lb 

REAL**  SIG(512) »  RNQIS ( 5 12 ) >3(512) * IMPULS ( 512 ) * SUMNOI ( 512 > 

LOG  I CAL*1  F ILNAM (15) *0TFLNM ( 15 ) *TEMP<6) , KEYPR * EXT ( 3 ) 

INTEGER  IX (2) 

DATA  IX/0.0/ 

CALL  JTITLE ( ' ADDNOI ' *6*  * 1 .2) 

TYPE  30 

FORMAT!/*'  Enter  time  data  input  filename  (used  to  compute  noise 
TYPE  35 

FORMAT ( '  statistics ) J  '**> 

CALL  GETDFN ( FILNAM ) 

IF ( OPNFIL  <  3 *  F ILNAM* ' R ' ) «NE ♦  ,TRUE.)G0T0  10 
CALL  GETDAT(3*NP*TFACTR*SIG* .TRUE.  ) 

DO  90  1  =  1*  NP 
S  <  I  >  - 1  * 

Compijt  statistics  for  input  file  data 


riiUS 


Input:  Reference  file 

Output :  Noise  only  file 


Input  J  Impulse  file 
Output:  Impulse  +  noise  Pile 


1 0  0 
100 


to 


*  ■.*  *■ 


CALL  TALLY ( SI G  *  3  *  TOTAL  *  AVERS  ?  3D3 1 G  *  VM I N  *  VMAX  ?  NP  *  1 > I£R? 


IF(IER.EG.0)G0T0  105 


TYPE  1 0 0  > I E R 

FORMAT!/* '  ERROR  NO.  '*12*'  IN  COMPUTING  STATISTICS 
GO  TO  200 

TYPE  110* (FILNAM( I ) *  1  =  1  *  15) 


FORMAT-://*'  STATISTICS  FOR  INPUT  FILE:  '*15Al*/> 


) 


TYPE  120* AVERS* SDSIG 

FORMAT < 3X *' Ave rage  =  '*012.3*'*  Standard  deviation  =  ' *  G 1 2 . 5 ) 
TYPE  130 

FORMAT! //*' -Enter  SNR  (relative  to  input  file  data)  to  compute 
noise  data  J  ' ) 

TYPE  131 

FORMAT ('  (INTEGER  <100)  '*$> 

ACCEPT  *  *  DBDUM 
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DBDWN=A3S( BBDWN ) 

3DN0JS*EXP'AL0G(10.0)*(AL0G10<3DSIG>-DBDUN/20.0> ) 

,* 

Generate  random  number  date  with  a  daussian  distribution 
C  with  standard  deviation*  SDNQI3*  and  mean  -  0.0  (assumes 

C  input  file  has  no  DC  offset 

C 

TYPE  *>  '  ' 

TYPE  -K»  '  Now  computing  random  generator  seed.' 

TYPE  **'  Depress  any  key  to  continue.  ' 

11  CALL  RANDU <IX(1)»IX(2>* RNO I S ( 1 >  > 

CALL  CHECK( IDUMMY > KEYPR > 

IF ( KEYPR  ,  HE .  .TRUE,)  GOTO  11 
DO  140  I-lrN'R 

CALL  GAUSS<IX*SBN0IS*0.0*RN0IS<I)  ) 

140  CONTINUE 

C  Get  statistics  for  noise  data 

CALL  TALLY <RN0IS»S» TOTAL *AVERN*SDNOIS*VMIN* UMAX *NPr 1  * IER) 

IF< IER.EG.0>G0T0145  ' 

TYPE  100* IER 
GO  TO  200 

145  3DR=20*AL0G10(3DN0IS/SDSIG) 

TYPE  150  '  > 

150  FORMAT ( // * '  NOISE  STATISTICS  t  '*/)  5 

SDR  ■  -SDR 

TYPE  1'60*AVERN>SHN0IS*SDR 

-40  FORMAT  <5X* 'Average  =  '*G12.5>'*  Standard  deviation  =  '»G12.5»/ 

1*5X* 'Actual  SNR  =  '»G12.5»'  dB'> 

EXT ( 1 )  =  '  N ' 

C  IDBDUN  -  INT(DBDWN)  +  100  !  needed  since  ITOA  pads  w/  null  char 

IDBDWN  =  INT(DBDWN)  !  Convert  SNR  to  text 

CALL  IT0A< IDBDWN*TEMP)  !  Get  SNR  in  dB's  into  TEMP- 

EXT(2)  =  TEMP ( 5 )  !  Put  SNR  in  dB's  at  last  of  filename 

EXT ( 3 )  =  TEMP <  6 )  !  first  4  chars,  are  stripped  off 

CALL  FIL EXT (FILNAM* EXT) 

TYPE  170*  FILNAM 

170  FORMAT  (/>  '  Attempting  to  write  noisy  dat3  to  '*15A> 

IF(0PNFIL(2*FILNAM* 'W' >  ♦  EG  *  .TRUE.)  GOTO  260 
230  TYPE  270 

270  FORMAT ( /  *  '  Enter  output  filename  t  '*$> 

CALL  GETDFN (FILNAM ) 

300  IF(0PNFIL(2*FILNAM» 'U ' >  .EG*  .FALSE.)  GOTO  280 

260  CALL  PUT DAT (2*NP*TFACTR*RN0IS) 

TYPE  #*'Data  written.' 

TYPE  #*  '  ' 

TYPE  210  ,  .  ,  , 

210  F0RMAT<//*'  Do  you  wish  to  add  noise  to  an  impulse  (system)  file 

1  '*«> 

IF  (ASK('N')  .EG.  .TRUE.)  GOTO  200 
225  TYPE  220 

220  FORMAT ( / * '  Enter  filename  for  impulse  data  t  '*$> 

CALL  GETDFN(FILNAM) 

IF  (0PNFIL< 3* FILNAM* 'R' )  .EG.  .FALSE.)  GOTO  225 
CALL  GETDAT  <  3* NP2  *  TFACTR  * IMPULS  * .TRUE. ) 

IF  ( NP2  .LE,  NP)  GOTO  227 
TYPE  228 

228  FORMAT ( / * '  Number  of  data  points  is  to  large!!  Try  again.') 

GOTO  225 
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•  HF'l 


SUMMOr-D  -  IMPULS'D 
0  0  '<  T 1 1 )  !J  £ 

*'  U-  PILE  X  r  •  FILM  A  M  ,  £  ;<  T  > 
i"'D  1  DD  F ILUAH 
■ '  ' 1  i !  A T  •'  /  f  '  A  t  l (J Hi f>  t  i  i  r iJ  to 
IF  <  CPNFIL  <  4  *  FILNAM  »  ‘ID  • 
CY  PE  222 


r  N  0  I  3  (  I 


u  r  lie  i  Hi  p  u  L  j  <j  n  o  i  <  y 
.EG.  .TRUE.)  GOTO  22:'. 


F ' 1 K  M  A  T  ( 1 X  «  '  £  n  t  *-  r  f  i  i  e  n  a  mo  J  '  >5) 
GALL  G E  T  D  F N  (  F I L  N  A M  > 

GOTO  2  X 

C  ALL  P  U  T  D  A  T  (  4  y  :•!  P2 > TFACTR » 3UMN0I  ; 

T  V  P  E  D  'Data  stored.' 

CONTINUE 
CALL  EXIT 

E  >!  D 


-  J  j  t  -3 


to 


SUBROUTINE  TALLY 
PURPOSE 

CALCULATE  TOTAL r  MEAN  >  STANDARD  DEVIATION j  MINIMUM j  MAXIMUM 
FOR  EACH  VARIABLE  IN  A  SET  (OR  A  SUBSET)  OF  OBSERVATIONS 

USAGE 

CALL  TALLY ( A  y  S  »  TOTAL y  AVER , SD  y  VM IN  y UMAX  y  NO  y  MV y I ER ' 

DESCRIPTION  OF  PARAMETERS 

A  -  OBSERVATION  MATRIX?  NO  BY  NV 

3  -  INPUT  VECTOR  INDICATING  SUBSET  OF  A.  ONLY  THOSE 

OBSERVATIONS  WITH  A  NON-ZERO  S(J>  ARE  CONSIDERED. 
VECTOR  LENGTH  IS  NO. 

TOTAL  -  OUTPUT  VECTOR  OF  TOTALS  OF  EACH  VARIABLE.  VECTOR 
LENGTH  IS  NV. 

AVER  -  OUTPUT  VECTOR  OF  AVERAGES  OF  EACH  VARIABLE.  VECTOR 
LENGTH  IS  NV. 

3D  -  OUTPUT  VECTOR  OF  STANDARD  DEVIATIONS  OF  EACH 
VARIABLE,  VECTOR  LENGTH  IS  NV. 

VM l N  -  OUTPUT  VECTOR  OF  MINIMA  OF  EACH  VARIABLE.  VECTOR 
LENGTH  IS  NV. 

V MAX  -  OUTPUT  VECTOR  OF  MAXIMA  OF  EACH  VARIABLE.  VECTOR 
LENGTH  IS  NV, 

MU  -  NUMBER  OF  OBSERVATIONS 

NV  -  NUMBER  OF  VARIABLES  FOR  EACH  OBSERVATION 

IER  -  ZERO-  IF  NO  ERROR. 

-  DIF  S  IS  NULL.  V  M I N  =  1  .  E  3  7  y  V  M  A  X  =  - 1  .£37.  y  S  D  =  A  V  E  R  =  0 

-  2y  IF  S  HAS  ONLY  ONE  NON-ZERO  ELEMENT.  VMIND'MAX. 

S  D  =  0 . 0 

REMARKS 

NONE 

SUBROUTINES  AND  FUNCTION  SUBPROGRAMS  REQUIRED 
NONE 


METHOD 

ALL  OBSERVATIONS  CORRESPONDING  TO  A  MON-ZERO  ELEMENT  IN  S 
VECTOR  ARE  ANALYZED  FOR  EACH  VARIABLE  IN  MATRIX  A. 
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W  -T-*  C<« 


c 

r* 


C 


TOTALS  ARE  ACCUMULATE!!  AND  MINIMUM  AND  MAXIMUM  VALUES  AK- 
f-OUND.  FOLLOWING  THIS,  MEANS  AND  STANDARD  DEVIATIONS  ARE 
CALCULATED.  H  E  DIVISOR  F  G  R  3  T  A  N  D  A  R  D  D  E  V I A  T  1 0  N  IS  ONE  ■  -  3 

THAN  THE  HUMBER  OF  OBSERVATIONS  USED. 


UUSROUTINE  TALLY ! A . S, TOTAL » AVER , SD , VMIN  > UMAX > MO ,NV ,  JER) 
DIMENSION  A  !  1  >  *  3  t  1  ) , TOTAL! 1 ) , AVER! 1 > ,SD! 1 ) ,  V  M I N  ( 1 > , UMAX! 1 > 

CLEAR  OUTPUT  VECTORS  AND  INITIALIZE  VMIN, UMAX 

IER-0 

DO  1  K - 1 r  M V 
TOTAL ( K ) =0 . 0 
AUER (K) =0.0 
3  D  \  K  >  =  0  ♦  0 
UM I N ( K  >  =  1 . 0E37 
UMAX  !!<)=-!  .0E37 

TEST  SUBSET  VECTOR 


3  C  N  T  =  0  *  0 
DC  ?  J=1,N0 
I J= J-NO 

I F ! S ( J  5 )  2,7,2 

2  SCNT--SCNT  +  1 .0 

r 

w 

C  CALCULATE  TOTAL,  MINIMA,  MAXIMA 

i  * 

DO  6  1=1, NV 
I J=I J+NO 
;<=A(IJ) 

TOTAL! I)=TOTAL(I)+X 
IF ! X-VMIN ! I ) )  3,4,4 

3  UM I N ! I ) =X 

4  IF ! X-VMAX! I )  )  6,6,5 

5  UMAX ! I ) =X 

6  SO ! I ) =SD  !  I ) tX*X 

7  CONTINUE 
C 

C  CALCULATE  MEANS  AND  STANDARD  DEVIATIONS 

/•% 

IF  I SCNT ) 3 , 8  >  9 
S  IER=1 
GO  TO  15 
9  DO  10  1=1 ,NV 

10  AVER ! I ) “TOTAL ( I ) /SCNT 
IF  ( SCNT-1 ♦ 0 >  13,11,13 

11  IER=2 

00  12  1=1, NV 

12  SD ( I) =0 ♦ 0 
GO  ro  15 
DO  14  1=1, NV 

SD!  I)=SQRT! ABS! ! SD < I ) -TOTAL ! I ) *TOTAL ! I ) /SCNT ) / ! SCNT-1 . 0  )  >  ) 
RETURN 
END 


i 


“*• .-C  i  ,  >,»  */.'**  *<„ 


c 

c 

r 


SUBROU  f  ih'E  -  GAUSS 
PUR POSt 

COMPUTES  A  NORM ALL"  DISTRIBUTED 
MEAN  AND  STANDARD  DEVIATION 


RANDOM  NUMBER  WITH 


USAGE 

CALL  G  A  U  S S  < i X  j  S  *  A M  ? V ) 


DESCRIPTION  OF  PARAMETERS 

IX  -IX  IS  AN  INTEGER  ARRAY  OF  LENGTH  2.  THE  INITIAL  ENTRIES 
IN  THE  IX  ARRAY  SHOULD  BE  ZERO.  THEREAFTER »  IT  'JILL 
CONTAIN  PART  OF  A  UNIFORMLY  DISTRIBUTED  INTEGER  RANDOM 
NUMBER  GENERATED  BY  THE  SUBROUTINE  FOR  USE  ON  THE 
NEXT  ENTRY  TO  THE  SUBROUTINE. 

3  -THE  DESIRED  STANDARD  DEVIATION  OF  THE  NORMAL 
DISTRIBUTION. 


AM  -THE  DESIRED  MEAN  OF  THE  NORMAL  DISTRIBUTION 
V  -THE  VALUE  OF  THE  COMPUTED  NORMAL  RANDOM  VARIABLE 


c  REMARKS 

■:  THIS  SUBROUTINE  USES  RANDU  WHICH  IS  MACHINE  SPECIFIC 

SUBROUTINES  AND  FUNCTION  SUBPROGRAMS  REQUIRED 
C  RANDU 

C 

C  METHOD 

C  USES  12  UNIFORM  RANDOM  NUMBERS  TO  COMPUTE  NORMAL  RANDOM 

C  NUMBERS  BY  CENTRAL  LIMIT  THEOREM.  THE  RESULT  IS  THEM 

C  ADJUSTED  TO  MATCH  THE  GIVEN  MEAN  AND  STANDARD  DEVIATION. 

C  THE  UNIFORM  RANDOM  NUMBERS  COMPUTED  WITHIN  THE  SUBROUTINE 

C  ARE  FOUND  BY  THE  ROWER  RESIDUE  METHOD. 

r 


SUBROUTINE  GAUSS ( I  X > 8 , AM r V > 
DIMENSION  I X  <  2 ) 

A  =  0.0 

DC  50  1=1,12  ' 

CALL  RANDU < IX < 1. »  t IX(2)»Y) 

50  A-A+Y 

V  (  A  -  6 . 0  )  TC  S  +  A  M 

RETURN 

END 
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£  MODIFIED  FOR  RT-11  V4  22-MAR-84 

PROGRAM  IMPULS 

J***  ******* *#** ******************************************************** 
C  THIS  PROGRAM  ALLOWS  ONE  TO  GENERATE  AN  IMPULSE  TRAIN  WITH  ARBITRAY 
C  NUMBER  OF  SPIKES. THE  CREATED  IMPULSE  DATA  IS  IN  TIME  DOMAIN. 

C****** ****************************************** ********************** 
DIMENSION  DATA ( 1 024 ) 

LOGICAL*!  FILNAM< 13) 7  IANS 
C 

10  CONTINUE 
0  TYPE  16 

16  FORMAT  ( 1 X  7  '  ENTER  IMPULSE  DATA  FILENAME:  '»*> 

CALL  GETDFN  <  FILNAM) 

IF  (0PNFIL(3»FILNAM» '  U '  )  .NE,  .TRUE.)  GOTO  5 


TYPE  30 

TO  FORMAT \ / »  '  SELECT  NUMBER  OF  DATA  POINTS  (2**1) t  '»£) 

ACCEPT  *  7  HP 
RNP=NP 
WRITE ( 3 ) RNP 
TYPE  40 

40  FORMAT ( / r '  HOW  MANY  SPIKES  IN  FULL  RECORD?  '»$) 

ACCEPT  *  7  NS 
TYPE  50 

50  FORMAT (/ > '  ENTER  DELTA  TIME  FACTOR  ( DEFAULT* . 80321 2S5E-S ) J  '?$) 
TMFCTR=.8032128514E-8 
ACCEPT  * 7 TMFCTR 
bJRI  TE  (  3  )  TMFCTR 
DO  <40  1  =  1 7  NP 

60  DATA ( I ) =0 . 0 

DO  100»I=l7NS 
TYPE  90 

90  FORMAT ( '  ENTER  SPIKE  LOCATION  (N)  AND  MAGNITUDE  (REAL)  I  '7$) 
ACCEPT  *  7  N  7  DATA  <  N ) 

100  CONTINUE 

DO  110  1  =  1 7  NP 
WRITE(3) DAT  A ( I ) 

110  CONTINUE 

200  CLOSE (UNI T  =  3  7  D ISPOSE= ' KEEP ' ) 

TYPE  210 

210  FORMAK/7'  WANT  TO  CREATE  ANOTHER  FILE?  (Y  OR  N)J  '7$) 

ACCEPT  2207IANS 

220  FORMAT ( A1 ) 

IF ( IANS  . EO .  ' Y ' ) GOTO  10 
END 
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I*-**-********-***-***'*****} 


PROGRAM  TRNSFM 

THIS  PROGRAM  PROVIDES  FOR  DISCRETE  FOURIER  TRANSFORMATION  OF 
DATA  IN  A  USER  INTERACTIVE  MODE. OUTPUT  DATA  IS  AVAILABLE  IN 
IN  AN  UNFORMATTED  FORM, THE  DATA  IS  HEADED  BY  NUMBER  OF 
POINTS  AND  SEPARAT  ION  BETWEEEN  FRQUENCY  PO I  NTS , FOLLOWED  BY  EXPON 
ENTIAL  WEIGHT  USED. THE  TRANSFORMED  DATA  IS  AVAILABLE  WITH  REAL 
PART  IN  ODD  RECORDS  AND  IMAGINARY  PART  IN  EVEN  NUMBERED  RECORDS. 


1100 


THIS  DATA  IS  IN  A  FORM  SUITABLE  FOR  PLOTTING  ON  THE  HP  PLOTTER  1 
7225A  USING  'FRQPLT'  SUBROUTINE.  | 

SUBROUTINES  USED  ARE  'FAST'  WHICH  IS  A  COMPLETE  PACKAGE  PROVIDING 
FOR  BOTH  FORWARD  AND  INVERSE  TRANSFORMATION  .OF  DATA. NOTE  HOWEVER!* 
THAT  IN  ITS  PRESENT  FORM  THE  PROGRAM  REQUIRES  RE AL * FUNCT I ON  OF  f 
SINGLE  VARIABLE)DATA  INPUT. THE  FFT  SIZE  ALLOWED  IS  1024  ALTHOUGH? 
THE  FAST  ROUTINE  PROVIDES  FOR  UPTO  4096  POINTS. FAST  IS  AVAILABLE 
AS  A  LIBRARY  SUBROUTINE  PACKAGE. 

DIMENSION  X  *  1 030 ) ,AUX*1030) 

L0GICAL*1  FILNAM  *18),  IANS 

DO  5  1  =  1  , 1030 
X(I)=0.0 
AUX ( I ) =0 . 0 
CONTINUE 
TYPE  15 

FORMAT*/,'  FREQUENCY  DATA  OR  TIME  DATA? * F  OR  T>',$> 

ACCEPT  1100, IANS 
FORMAT ( A1 > 

IF ( IANS  , EQ . ' F ' ) GOTO  120 
IF( IANS  .NE. 'T' )  GO  TO  16 
CF( IANS  .EQ. 'T' >  GO  TO  17 

TYPE  *,'  WRONG  DATA  TYPE, TRY  AGAIN?  ' 

GO  TO  10 


IF( IANS 
[F( IANS 


’  F ' ) GOTO 
’T'>  GO 
T '  >  GO 


OTO  120 
GO  TO  16 
GO  TO  17 

WRONG  DATA  TYPE, TRY 


AGAIN? 


CONTINUE 


OPEN  TIME  DATA  FILE 


1200 


1300 


TYPE  20 

FORMAT*/,'  ENTER  INPUT  TIME  DATA  FILE  NAME*.  ',*) 
CALL  GETDFN  *  FILNAM ) 

IF  (0PNFIL(3, FILNAM, 'R ' )  .NE.  .TRUE.)  GOTO  13 

READ ( 3 ) RNP 

NP=RNP 

TYPE  1200, NP 

FORMAT*/,'  NUMBER  OF  DATA  POINTS  =  ',14,/) 

READ ( 3 ) XFCTR 
TYPE  1300, XFCTR 

FORMAT*/,'  DATA  SPACING  IS  J',G12.6) 

DO  SO  1=1, NP 
READ*3,ERR=55)X(I) 

CONTINUE 
GOTO  60 

TYPE  READ  ERROR,  FILE  SCREWED  UP' 
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A  v- 

C 


C 

r\ 

U 

c 

70 


74 

C 


Q0  TO  200 

CLOSE ( UNI T  =  3, DISPOSE* 'KEEP ' ) 

PROMPT  THE  USER  FOR  EXPONENTIAL  WEIGHTING  FACTOR  WHICH  SHOULD  BE 
AS  CLOSE  AS  POSSIBLE  TO  1,0  <0.9999  IS  A  GOOD  CHOICE  FOP  Nnisc 
FREE  DATA). 

EW=1 .0 
TYPE  AS 

FORMAT  <  / »  '  ENTER  WEIGHTING  FACTOR  TO  BE  USED  <  <1,0  >:',*) 

ACCEPT  *,EW 

CALL  EXPWAT<X,EW»NP, 1 ) 

DEFINE  FFT  SIZE  AND  CHECK  FOR  CORRECT  NFFT  TO  AVOID  ALIASING. 
TYPE  72 

FORMAT ( / »  '  ENTER  DFT  SIZE,  LENGTH  MUST  BE  A  POWER  OF  TWO  J'S) 
ACCEPT  *  > NFFT 

IF'NFFT.LT.NP)  GO  TO  73 
IF(NFFT.GE.NP) . GO  TO  74 

TYPE  DFT  LENGTH  MUST  BE  >  NUMBER  OF  DATA  POINTS:* 

GOTO  70 

CONTINUE 


FMAX=1./<2.#XFCTR> 

FFCTR=2.*FMAX/FL0AT(NFFT> 

GET  READY  TO  COMPUTE  FFT  USING  "FAST" 


FORM  SEQUENCES  FROM  X(N)  AND  N*X'N>  FOR  DFT. 


90 

n 

w 

C 


c 

c 

C 


f* 

*4 


c 

c 

104 

105 


no 


DO  90  I  =  1 »  N  P 
AUX<I)=FL0AT<I-1)*X<I> 

CGNTINUE 

COMPUTE  FFT 

CALL  FAST ( X » NFFT ) 

CALL  FAST ( AUX,NFFT ) 

STORE  THIS  FREQUENCY  DOMAIN  DATA,  | 

OPEN  A  NEW  FILE  WITH  GIVEN  NAME, THE  FIRST  RECORD  CONTAINS  NUMBER  0 
FREQUENCY  POINTS  AND  FREQUENCY  SPACING. SUBSEQUENT  RECORDS  CONTAIN 
REAL  AND  IMAGINARY  PART  OF  TRANSFORMED  DATA , FOLLOWED  BY  REAL 
SIMAGINARY  PART  OF  N  TIMES  TRANSFORMED  DATA. 

TYPE  105 

FORMAT*/, '  ENTER  OUTPUT  FILE  NAME < FREQUENCY  DATA)'?*) 

CALL  GETDFN ( FILNAM ) 

IF  <  OPNFIL  <  2.,  FI  LMAM ,  'W'  )  .  NE  .  .TRUE.)  GOTO  104 
R£AL=l.+FL0AT<NFFT)/2. 

WRITE<2)REAL»FFCTR»EW»FMAX 

DO  110  1  =  1, NFFT  +  1,2 

WRITE(2)X<I) , X  < I  + 1 ) ,AUX<I ) ,AUX<I  +  1) 

CONTINUE 

CLOSE* UN IT=2, DISPOSE*' KEEP' ) 
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tff  ■£> 


!;U  TO  200 


PROCESS  1 NG  FREQUENCY  DOMAIN  DATA 
TYPE  122 

FORMAT  .  /  »  '  ENTER  FREQUENCY  I?  AT  ft  FILE  NAME  •  INPUT  KATA)  .*•,?> 
CALL  GETDFN ; FILNAM  > 

IF  ? OPNFIL ( 3  »  FILNAM  »  '  R '  >  . NE.  .TRUE.)  GOTO  120 

READ ( 3 > RNP , FF CTR »  EU »  FMA  X 

NP-RNP 

IF (EU.EQ.O.O)EU-l .0 
TYPE  1 200 , NP 
TYPE  1300 » FFCTR 
TYPE  1400rEU 

MOO  FORMAT ( / »  '  EXPONENTIAL  WEIGHT  USED  IS  :',F6.4) 

DO  140  1=1 »NP 

READ ( 3 »  ERR=55 )  X  {  J )  , X ( J+ 1 >  » AUX C J >  ,  AUX?  J+l  ) 

140  CONTINUE 

CLOSE <  UN  I  T  =  3,  DISPOSED  KEEP' ) 

150  TYPE  72 
ACCEPT  *»NFFT 

IF?NFFT.LT.2*<NP-1> )  GO  TO  151 
IF?NFFT.GE.2*<NP-1 ) )  GO  TO  152 

151  TYPE  *>'  DFT  SIZE  MUST  BE  >2*<NUMBER  OF  FREG  POINTS-!)' 
GOTO  150 

152  CONTINUE 
TMAX=1 ./FFCTR 


PERFORM  INVERSE  TRANSFORMATION  TO  RECOVER  REAL  DATA  FROM 
FFT  TRANSFORMED  DATA. 


170  CALL  F38T ( X »  NFFT ) 

CALL  EXPWAT(X»EW»NFFT,1) 

C 

C  STORE  THIS  DATA 

C  OPEN  A  NEW  FILE  WITH  GIVEN  NAME. THE  FIRST  RECORD  CONTAINS  NUMBER 

2  DATA  POINTS  AND  SECOND  CONTAINS  SAMPLING  INTERVAL  FOLLOWED  BY  TH 

C  DATA. USE  ' HFPLQT '  TO  PLOT  THIS  DATA. 

C 

TYPE  175 

FORMAT ?  / »  '  ENTER  OUTPUT  TIME  DATA  FILE  NAME  t',$> 

CALL  GETDFN  < FILNAM  > 

IF  (OP NFIL?3, FILNAM, 'W' )  . NE .  .TRUE.)  GOTO  174 

REAL=FLOAT<NFFT) 

XFCTR=TMAX /FLOAT (NFFT) 

WRITE? 2) REAL 
WR I TE  ?  2 ) XFCTR 
DO  130  1  =  1, NFFT 
WRITE (2)X(  I  ) 

CONTINUE 

CLOSE  (  UNI  T*=2  »  DISPOSE^  '  KEEP  '  ) 

2 'j 0  TYPE  .210 

2-.0  FORMAT?/,'  DO  YOU  WANT  TO  TRY  AGAIN?  <  Y  OR  N>:  '5? 

I A  N  3  =  '  Y  ' 

ACCEPT  1100, IANS 

[F ( IANS  .EG.  '  Y ' ) GOTO  1 


SUBROUTINE  EXPUAT ( X , A , N , K ) 
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DIMENSION  X(2> 

IF  ( A »  EO  # 1 . 0 ) RETURN 
FX=*1  . 

DO  10  £  J 1 »  N 

xm=x<i>*Fx 

:  F  <  K  .  i3T  .  0  >  QO  TO  13 
IF  < K ♦ LE#  0  5  GO  TO  20 
13  FX=FX*A 

20  FXsFX/A 

10  CONTINUE 

RETURN 
END 


I 


■  I'UGRAM  !-i I  ENER 


in  la  ncfj 


• e  l  -  o  w1 


HEEDS  - 
2  1 -Mur-85 


1 2  - '  f  A  Y  -  3  5 


Chuck  NFFT  <53*3 most  lard-eel  of  ( or- 1  *  .  ir-2  :• 
add  time  dome  in  plot  aft  a  r  Ira  n  s  f  u  r 
Changed  at  a r tind  p  t .  in  loop*  f  r o m  2  to  2 
and  chunked  X<  >  to  REF  (  )  f  Y()  to  FLAUO, 

I  to  REAL*  3  I M  A  G  to  make  program  Plow  u  a  a  1 u  r 
to  follow. 

Add  if  one  wants  to  plot  the  time  domain  result 
to  the  plotter. 


lo-Jul-35 


Corrected  YPLOMF  3  General  update 


THIS  PROGRAM  PERFORMS  WIENER  FILTERING  ON  THE  DATA  IN  FREQUENCY 
DOMAIN  ACCORDING  TO  THE  EQUATION  R ' < f ) *F < f ) / < R < F ) **2+FACTR ) , 
FACTR  IS  THE  DESENSITIZING  COEFFICIENT  AND  IS  USUALLY  SET  TO 
REDUCE  THE  ABOVE  DIVISION  TO  A  SMALL  VALUE  OUTSIDE  THE  RANGE  OF 
INTEREST, I ,E. , OUTSIDE  THE  FREQUENCY  BAND  OF  THE  TRANSDUCER. 


LOGICAL*!  FILNAM(IS) » REFFIL ( 1 5 ) ?  FLAUFL (  15 )  , ANSWER 
REALJK4  FFCTRt  XFCTR,  XFCTR2*RMAX,  CUTOFF  » PERCNT , TMFCTR 
I M T  E  G  E  R  #  2  REAL  ,  I  MAG »  I »  NF'  1  ,  NP2  ?  NP3 
DIMENSION  REFt 1030) t FLAU< 1030) 

CALL  jT I TLE( 'WIENER' ,0, ' I ',2.1) 


DO  20  1=1,1030 
REF ( I )=0. 0 
FLAW ( I >  =0 . 0 
CONTINUE 

SET  DATA  FROM  FILES 

CALL  GETD ( REFFI L , FLAUFL  r  NP1 »  NP2  »  XFCTR  »  XFCTR2 , REF , FLAW ) 
r  M  F  C  T  R  =  X  F  C  T  R 


PROMPT  THE  USER  FOR  EXPONENTIAL  WEIGHTING  FACTOR  WHICH  SHOULD  3: 
AS  CLOSE  AS  POSSIBLE  TO  1.0  (0.????  IS  A  GOOD  CHOICE  FOR  NOISE 
FREE  DATA). 

EU-i . 0 
T Y PE  30 

FORMAT-;/,  '  Enter  weighting  factor  to  oe  used  <  11,0  >t  *5. 

A  C  C  E  P  T  *  >  E  W 

CALL  EXPWAT ( REF , EW , NP1 , 1 ) 

CALL  EXPWAT ( FLAW , EW  t NP2  *  1 ) 

DEFINE  FFT  SIZE  AND  CHECK  FOR  CORRECT  NFFT  TO  AVOID  ALIASING. 
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F F  C T R - 2  ,  *  F M  A  X / FLOAT  (  M FFT  > 
COMPUTE  FFT 


!  i  r  _  >■  ? 

T  V  P E  $  >  '  CrsnsPoriTii  n 3  R <5  f  e  9 n ce  d a  1 o  t  o  3 - ?  1  a n e  ♦ 


CALL  r  A  3  T  (REF  ?  N  F  F 7' ) 

T Y PE  *  ?  '  T  rens  fo  rising  Fla w  d a t a 

\i.-  t  t 

+•  f 

FAST < FLAW >NFFT) 


rU  0“lififii 


•  i  t”  IE. 

C  A  L  L 


r  0  R  M  A  T  (  /  ? 

b  u  s  u  U 

wish 

to  plot  the 

IF  CASK-'N 
TYPE  90 

’  )  .EG. 

T  -•  t « c 

.  >  GOTO  120 

FORMAT- /»  ' 

D  c  ■=•  o  u 

wish 

to  Plot  the  i 

ACCEPT  100 
F U  R M AT ( 1 A ) 

» ANSWER 

IF  (  ,  N u  T  .  C 

A  N  S  W  E  R  * 

EG  ♦  '  F 

. OR . ANSWER . E 

IF  ( ANSWER 

»EQ. ' F  ' 

)  CALL 

YPLCMF { FLAW 

IF  'ANSWER 
TYRE  110 

.  SQ .  '  R  ' 

)  CALL 

YPLOMFCREr > 

FORMAT  < / r  ' 

A  ri  o  t  h  e 

r  plot 

'  7  %  ) 

IF  '  A  3  K  ( ' N 
CONTINUE 

'5  .HE. 

.  TRUE 

, )  GOTO  30 

r  e  a  »  c  c  m  a  i  ri  a 


o r  ri, 


WIENER  TRANSFORM 


-LAN •  1 )  =0 . 0 

•  A  I  I  O  i  « 

I*  L  H  t*)  >  ril  /  “  V  *  ‘J 

REF  •  1  ;  -  0 . 0 
REF ( 2) =0 . 0 


:  0  v  u  ♦  C  ♦  C  O  :Ti  P  O  fi 


30  REAL-3  »  MFFT  *  2 


1  M  A  3  -  R  L  A  L  +  1 

E  P  1=RE  F  •  R  E  A  L  ;  *  FLAW  <  READ +REF  •  IM  A  G  )  H  F  L  A  W  <  I  M  A  3 
T  E  M  P  2  =  R  E  F  •  REAL  )  *FLAW  <  IMAG )  -REF  •  I  MAG  >  *FLAN  i  REAL 
FLAW(  READ-TEMPI 
FLAW'  I M  A  G  )  =  T  E  M  P  2 
CONTINUE 


COMPUTE  SQUARED  MAGNITUDE  OF  REFERENCE  DATA 


j"i  x  —  r» 


DO  140  REA 


r  \/  *i  p 

i  h  U 


-  FV  22  H  L- 


3  7  N  FFT  7 
P  3. 


•  R E A L  )  - R E F  :  R E  A L  *  * 2  +  RE-  IM A 3  •  P * 2 


MAY  .  L 


iL  )  )  R h  A  X  =  R F  -  !-: E  *  L 


lL  Lli 


st'ifiy  r-3 p i- :"i  L 3 b  .  V ersio n  1.0  3  ! 

L  : SICA L  S'  1  M C H A R  ,HCH A R 
DIMENSION'  X  <  3 1 5  )  »  Y  (  5 1 5  )  >  Y F  I  T  (  5  1 5  ) 

REAL  A<  50 »  50 ) »  B ( 50 ) »  Z ( 200  >  »  C ( 50  >  »R(51) » P (500 ) - D  <  5? 
LOGICAL *1  IAN'S 
INTEGERS  TEKCRT  * CRTAGN 
DATA  TEKCRT/ 3/  CRTAGN/4/ 

CALL  JTITLE-;  '  TEKDEC'  ?6,  '  =  '  ,  1.0) 

CALL  CLEAR 

GET  INPUT  DATA  FILE 


DO  ?  1=1 >256 
YF  IT! I ) =0 . 0 
Y ( I ) =0 . 0 
CONTINUE 

CALL  INPUT  <  Y  >  N P  >  T F C T R  > 1 ) 

IF  ( NF  . LE.  256)  GOTO  15 

TYPE  *  >  Number  of  data  points  must  be  <  =  256 
GOTO  16 

CREATE  TIME  AXIS 


DO  20  1=1 fNP 
X ( I > = FLOAT ( 1-1 ) 

CONTINUE 

PROCEED  TO  GET  PARAMETERS  FOR  SPLINE  FITTING. 


TYPE  40 

FORMAT! />'  E  n  t  e  r  ratio  of  k  n  o  t - 3  p a  c i n S  to  d a  t  a  s p  . 
ACCEPT  * >  K R 
00  50  1=1 >200 

2  (  I  >  =FLOAT  (  I  - 1  )  AFLOAT  ( !<R  ) 

CONTINUE 
TYPE  60 

F 0 F: MAT  (  /  >  Enter  order  of  splines)  ‘  ?  %  ) 

ACCEPT  % f K  K 
TYPE  70 

F 0 R M A T  (  /  f  '  E 11  ter  n umber  of  basic  splines)  '  •  5  ) 
ACCEPT  *?N 
KK  1  =KK-  1 

00  ?0  1  =  0?  K K 1 
/•»  .  ' 

-.3  -  U  * 

LU=KR*(KK-I ) 

DO  30  L  =  0  >  L U 
T 1  =  F  L  0  A  T  (  L  f  K  R  *  I  > 

T  2  =  F  L  0  A  T  (  L  ) 

3=  3  + BSP < T1 >  Z  ? 1 >  KK ) * 3SP( T2  >  Z  > 1 . KK ) 

CONTINUE 

•• :  I  +  1  )  =  S  : 

CONTINUE 


r  ..  • « "i  >  »■»  j 


CONTINUE 

no  140  J  =  1 ?  N 
e  •:  j  >  =  o . 

M L  =  K R #  <  J-l )  +  l 
fi  U  -  K  R  'fi  J-1t  K  !\  )  + 1 
00  130  I=ML?MU 
JKL 1  =  J 

B ( J ) =  B  < J ) +Y ( I ) *BSP ( X  ! I > 7 Z • JKL 1 > KK ) 
CONTINUE 
CONTINUE 

!\3=KK-1 

CALL  BW3 ( A » N  ♦ KB »  3 » C ) 

DO  150  1  =  1  »  NP 
T  =  F  L  Q  A  T  ( I  -  1 ) 

Y  F  I T  (  I  )  =  P  3  <  7  ?  Z  «  C  ?  N  •  K  K  > 

CONTINUE 
DO  155  1=1 ,NP 
X  ! I )  =  F  L  0  A  T ( I - 1 )  *  T  F  C  T  R 
CONTINUE 
TYPE  160 

F 0 R M A T ( / r  '  Wish  to  plot  Tit  to  reference  date? 
IF  (  A  3  K  <  '  Y  '  )  ♦  NE  .  ♦TRUE.)  GOTO  11 
CALL  CLEAR 
TYPE  ' 

TYPE  Plotting  reference  data...' 

CALL  PLOTXY!XrY?NP? TEKCRT) 

TYPE  't  >  '  Plotting  spline  fitted  data...' 

CALL  PLOTXY ( X  *  Y  F I T  r  N  P » CRT  A G  N ) 

TYPE  ISO 

FORMAT!/?  '  Do  sou  wish  to  r  e  e  I  o  t  this  si  r  a  p  h  ?  ' 
IF  (ASK! ' N ' )  ,EQ.  .FALSE.)  GOTO  170 

BEGIN  DECONVOLUTION  STEPS. 

DO  190  1  =  1 ?  N P 1 

Y  <  I )  =  0  . 

GET  SCATTERED  DATA  FILE  NAME 

CALL  INPUT! Y»NP1 ?TFCTR?2) 

IF  ! N  P  . LE.  256)  SOTO  196 

TYPE  Hr'  Number  of  data  points  must  be  -1=  25- 
GOTO  195 

GET  PARAMETERS  FOR  SCATTERED  DATA  FITTING. 

T  'PE  210 


F  0  R  M  A  T  !  /  i 

E  n  t  e  r 

ACCEPT  * ?  K 

n 

TYPE  220 

FORMAT ( / ? ' 

E  n  t  e  r 

ACCEPT  t 7 L 

KN=KK+K3 

L  L  =  1  M  P  K  N  - 1 

>  %  K  R  - 1 

'  r  !  1  (:T 

-  r  ,  •  j  .  ♦  _ 

0  0  >  T  Y  P 

3  --  b  t  C  (  L  i 
CON  TIN'JE 


CONTINUE 
L  L  N  r  i\  S  -  2 

IF  (  KRX  (  N-rL  +  KM-I )  .  GE  .  H  )  TYPE*  »  '  (  NtLtKN-1  )  K K R 
DO  260  1=0  ? LL 
3  =  0  . 

LI=(N+KN-I-1 )*KR-1 
DO  250  LS  =  1 ?  L I 
3»S+P ( LS ) *P ( L3+ • *KR ) 

CONTINUE 
R  < I  +  1  ) =  S 
CONTINUE 
DO  290  1  =  1  ,  L 
DO  290  J  =  I ?  L 
Hi  If J) =0,0 
IF (LL.GE.L-1 )LL=L-1 
DO  310  I =  0  j  L L 
L I  =  L- 1 

DO  300  J  =  1 ?  L  I 
A<  J,  J+I  >=R(  1  +  1  ;• 

JI=J  +  I 
CONTINUE 
CONTINUE 

KL=<N+KN-1 )*KR-1 
DD  350  K 1  =  1 ?  L 
S-0  > 

K  L  1  =  K  L 

KL2=KL+(K1-1 >*KR+1 
IF <KL2.GE.NP1)KL1=NP1-(K1-1 >*KR-1 
DC  340  1=1? KL1 
S=S+Y( If (Kl-l )*KR+1 >*P< I > 

CONTINUE 
3  (  K 1  )  =S 
CONTINUE 
IBW=N+KN-2 

r-  :  ISU  .  GE  .  L-l )  ISy  =  L-.t 
CALL  BUS ( A  ?  L  ? I BW  ?  B  ?  D  > 

DO  360  1  =  1  ?  N p 1 
T --.-FLOAT  -  1-1  ) 

Y  F  I  T  (  I  )  =  P  S  (  T  ?  Z  ?  D  ?  L  ?  k  3  ) 

TYPE  365 

FORMAT!/,'  Plot  computed  impulse  response? 
IF  (ASK('Y')  .  NE  .  .TRUE.)  GOTO  12 

CALL  CLEAR 

0 A  LL  PLOT  X  Y { X  ? Y  F I T  ?  N P 1 ?  T  E K  CRT) 
r  y  p  :r  i  a  r. 

!.'  F  •:  ASK  •  '  N  '  )  .  ME  .  .  TRUE  .  )  GOTO  3?Q 

TYPE  37 5 

F 0 R MAT ( / • '  Wish  to  store  impulse  re s ponse? 
IF  <  ASK ( ' Y ' >  .NE.  .TRUE.)  GOTO  390 
TYPE  330 

F  0  P  M  A  T  '  ■'  ?  '  Enter  i  m  '■  u  1  s  e  r  e  s  p  o  n  s  e  d  s  t  o  Til 


■1 1  $ 
4  o 


0>  59  0 

TYPE  600 

$  4  C  0 

FORMAT ( / > 

IF  (ASK -I7 

i  «... 

CALL  EXIT 

R 

o  ;  Ur 

E  A  L! 

BTAST  COMPUTE!' I  OHS  OF  FIT  T0  FLAW  DATA. 

DC  no  I-l tMPI 
!_ L  -  (  N t K ii  —  1  >  * K R  -  1 

DO  4  00  J- 1  ♦  L 

IF ( I -1- <  J- 1  > *KR . LT . 1 >  GO  TO  400 
IF  <  I  -1  -  <  J-l  >  *KR  .  GT . LL  >  GO  TO  400 
S  =S  +  Ei  (  J  )  #P  <  I  - 1  -  0  J - 1  )  K K  ) 

CONTINUE 
VEIT • I ) =S 
TYPE  420 

FORMAT <  / »  '  riot  fit  to  flaw  data  ?  ' »  $ ) 

IF  (ASK('Y')  * NE.  .TRUE.)  GOTO  590 
CALL  CLEAR 

CALL  PLQTXY • XrY> NP1 > TEKCRT > 

CALL  PLQTXY  <  X  »  YF I T t  NP 1 » CRT AGN ) 

TYPE  130 

IF  < A S K ( ' N ' )  .ME.  .TRUE.)  GOTO  570 

Wish  to  p  e  r  f  o  r  to  another  miracle? 


$  ) 


y  «■ 


o  J  i?  r»uu  \  i  c.  u  u  i  o  z.  \  u  i  y  u  *-  7  iOui.  u  > 

Wridht-Patterson  Air  Force  Base 

The  subroutines  used  i n  this  proarsm  were  t a k e n  from 

-programs  for  Digital  Signal  Processing' 

IEEE  Press  *  197? 

345  East  47  Strett»  New  York*  NY  10017 
Sponsored  by  the  IEEE  Acoustics*  Speech  *  and 

Signal  P r  o  c  e  s  s i n d  Society 
Lib.  of  Co ri dress  Cat.  Card  *  79-89023 
IEEE  Book  *  0-37942-128-2  (paperback  ver.) 

t  0-37942-127-4  (hardback) 

Also  published  by  John  Wiley  %  Sons*  Inc. 

Wiley  Order  #  0-471-05961-7  (paperbac k  ver. ) 

I  0-471-05962-5  (hardback) 

DIMENSION  U  1(1. 024)  >1)2(1024) 

TW0PI=8.0*ATAN( 1.0) 

CALCULATE  A. UNIFORM  DISTRIBUTION.  SEED  FDR  U2  IS  UK 


CALL  UNIDSTPJl  * 1 024  *  I  SEED  ) 

I3EED=U1 ( 1024) *16334 
CALL  UNI D ST (  U  2  ?  1024  * ISEED) 

MOW  COMPUTE  A  NORMAL  DISTRIBUTION  FROM  U1  AND  U2. 
PLACE  RESULT  IN  Ul. 

DO  10  1=1 * 1024 

Ul(I) =SQRT(-2.0*AL0G(U1 ( I ) ) ) *C03 ( TW0PI*U2( I ) ) 
CONTINUE 
RETURN 
E  N  D 


SUBROUTINE  UN I DST ( U * N * I  SEED ) 
DIMENSION  U ( 2 ) 


IF  (ISEED.EQ.O)GOTO  20 
1 3  -=  1 3  E  E  D 
DO  10  1=1 7 N 

IS IS=M0D( 131*13 * 16334) 
U  (  I  ) =  F L 0 A T ( IS )/l 6334.0 
CONTINUE 


£r>J 


PROGRAM  CEPSTR 


THIS  PROGRAM  READS  AN  INPUT  DATA  FILE  FROM  A  DISC  PREVIOUSLY  DEFINE 
THE  DATA  CONSISTS  OF  CONVOLVED  SIGNAL  IN  TIME  DOMAIN . EXPONENTIAL 
WEIGHTING  IS  USED  TO  ELIMINATE  SPECTRAL  ZEROES . COMPLEX  CEPSTRUM  IS 
COMPUTED  USING  A  SET  OF  SUBROUTINES  ADAPTED  FROM  LITERATURE. 

PHASE  UNWRAPPING  IS  ACCOMPLISHED  THROUGH  THE  TECHNIQUE  DEVISED  BY 
TRIBOLST. INTERACTIVE  PLOTTING  IS  POSSIBLE  ON  H-P  PLQTERS  THROUGH  US, 
OF  THE  ' H F I S P P '  PLOTTING  PACKAGE. 


V e r s ion  1.1 
Version  1,2 


10- Jan-35 

24- Ju 1-35  Modified  PLOT/ 


EDIT  history: 

Copied  from  H P C E P 3  and  modified  to  its  present  form  Feb.  13*1? 

-  GWL 


COMMON  / DATA// <515 ) *CY< 515) 

COMMON  PI rTWOPI *THLINC» THLCCNfNFFTfNPTS* N*L*H* HI  *  D  V  T  M  N  2 
DIMENSION  A U X ( 5 1 5 ) 

LOGICAL *1  I  A  N  S • N  P  A  G  E  <  2  ) *IPZER(10) *  ROTATE 
LOGICAL*!  FILNAMC15) 

INTEGER  I T P R L  <  6 ) *  G W 
LOGICAL  ISSUC 

CALL  JTITLEC 'CEPSTR' »6» »1.2> 

INITIALIZE  BUFFERS  AND  CONSTANTS 


ROTATE  -  .FALSE. 

THLINC-1 .5  !  set  compare  values  for  p h a u n w 

THLCON- » 5  ! 

P  I  -  4 , 0  *  A  T  A  N  (  i  «  0  > 
fwOPI =2 . *P I 

GET  INPUT  TIME  DOMAIN  DATA. 


»•  r.  r~  c* 

.  i  r  L  j 

F 0 P MAT''/,'  Enter  input  t  i  ® e-domain  data  f  i  I  e n a m e  <  c o n v o  1  v e d 
lata)  *  r  $  ; 

CALL  GETDFN ( F I LNAM  > 

IF  ( 0  P  N  F I L ( 3  *  F I L  N A M  >  '  R  '  >  . NE .  .TRUE.)  GOTO  7 

CALL  GETDAT  (  3  *  NF’  *  TFCTR  *  Y  *  .TRUE.  ) 

F- CTR  =  1 ,/vTFCT REFLOAT ( N P  )  ) 


'HE  SIZE  OF  DFT  AND  CHECK  IT! 


STH  F 


r  R  0  F' 


n  m  £•  i  s 


28 


NF  FT  ■■  if  ROM:  <  512  ) 

C !': e c K  Tor  N FFT  too  1  3  r a  =3  . 

Irinr'FT  .GT.  512)  GOTO  30 

Check  I' or-  NFFT  not  a  f-ouer  of  two 
L0G2NP-0 
I  TEMP=NFFT 

I F  <  I  T  £  M  P  ,LE.  1)  GOTO  503 
I  TEMP- I  TEMP / 2 
L0G2NP=LQG2NP+1 
GO  TO  34 

MPP0U2  =2*#LQG2NP 
IF (NFFT  . NE .  NPP0W2  >  GOTO  30 
NFFT  must  be  3  power  of  two 

-:nfft=fi.oa  r  <  nfft  :■ 

00  35  I  ~  1  ?  N F F T -r 2  !  USE  AUX  FOR  PLOT  PURPOSES 

A U X ( I )  =  ( I - 1 ) # T F C T R  !  AUX  is  used  as  X  date 

CONTINUE 

r  m a  I i 0  e  da t a  upon  re  a u e  s  t 
TYPE  600 

FORMAT  (  /  >  '  Do  sou  wish  to  normalize  the  -lets  '  >$) 

I F  (  A  S.K  (  '  N  '  )  ,EQ.  .TRUE.)  GOTO  2 
CALL  NORMAL ( Y»NP • 1 ) 

P R G  MPT  THE  U  3  E  R  F  0  R  E  X  P 0  N E  N T I A  L  WEIGHTING  F  A  C  T  0  R , 
WEIGHTING  FACTOR  SHOULD  BE  AS  CLOSE  AS  POSSIBLE  TO  1.0. 


E U  =  0  •  ???? 

TYPE  36 

FORMAT  (/»  '  Enter  w  e  i  d  h  1 1  r.  d  factor  (  Default  =  <  9  0  9  ?  )  '  ) 

TYPE  640 

F  0  R  M  A  T  1  X  t  '  C  e  p  s  t  r  3 1  processing  time  is  decreased  wit  h  d  e  c  r  e  a  s  e  d 


F  0  K  M  A  . 

(  1  X  > 

s  i  3  h  t  i 

na  fa 

ACCEPT 

*  >  EW 

IF  (EW 

.EQ. 

IF  (  E L 

.LT. 

IF  (EW  ,EQ.  1.0)  GOTO  512 
IF  (EL  .LT.  1,0)  GOTO  39 
r  Y pc  i  5 

FORMAT'./? '  Exponential  weidhtind  is 


-  Y  =  1  , 

DU  37  I=1?NP 
Y (  I  )  =  Y ( I ) *FX 
FA -F X M E  W 
CONTINUE 


PLOT  HEIGHT 


c  r  r:  r  c  n  ri  a  t  . 


-  0 K M A T  (  /  ?  Do  you  w.s h  1 0  - i c t  the  w e  1  d h  t a d  d a  t a  ’  .»  5  ;■ 
IF  (ASK('N')  .EG.  ,  TRUE.)  GOTO  53 
CALL  CLEAR 

C A L L  PLO T X  Y  •  A U X  ?  Y  ?  M P  -  I ;  RT-11  P L 0  T 
•  Y P E  111 

FQ:-"1m  T  '  /  -  to  you  wish  to  save  weishteC  late  '  * 


CO  MI'S  =  -SOUNDS  -  60  ♦  O^FLOhT  vMINITS)  -  3600 , 0*F  LGAT  <  ] 
' '  P  E  6  *3  0  -  L-:  a  -  I H  0  u  R  S  *  M  E  N I  T  S  ?  3  C  0  N  D  3 

UK  MAT  <  /  •»  '  .errors!  frocsssins  of  the  visidh  ted  Us  to 
o  f  '  *  R  ' .  3  ?  .•'  r  '  reoui red  '»I2r'  hro  .  •  '  *  1 2  r  '  a  ins  -  .« 


r  .  •.  iv  /■  •%  .•>  . 

-r  J  .*  !  -•  •„  J  _> 


*  i  t  fit  1  n  S  r  .« 


r  E  3  3  U  C  )  0  U  T  0  5  0 
Y £  4  6 

U  H  A  f  (  ■’  ’  '  ?  h  333  U  f:  W  r  a  ?•  P  1  fi  3  f  3  1  1  e  d  ♦  Try  tj  S  1  H  O  S  i  O  W  S  I 
e  i  si  h  tinsi .  ’  > 

0  TO  2 

HASH  UNWRAPPING  SUCCESSFUL.  WRITE  OUT  POLARITY  tPHi- 
NO  EXPONENTIAL  WEIGHTING.  USED . 

YPE  51 » ISNXf ISFX»EW 

U  R  fi  A  7  <  /  ?  '  3  i  «s  n  =  »  1 2  » /  r  '  Lines  r  j?  h  3  s  e  -  ’  1 

/  7  E;:r-'oneri  t  i  si  weighting  used  -  '  >  £10 . 3  > / > 

0 o hi ute  to  tsl  energy  in  csrst" u m 

0  FENG -0  . 

0  60  I - 1 >  N  F  F  T 
T 0 T E N G  -  T G T E N 3 -r C Y  '  I  >  *CY  (  I  > 

A U X  (  I  )  - C  Y  <  I  >  s  3  v e  C  Y  ds  t 3  f  c  •  •  ' e n e s  l e d  u o 

ON T INUE 

Y  P  E  &  1  1  T  C  T  E  N  G 

JR-V-tT  ‘Si'  Total  energy  m  oe^s  t  rum  -  «  £12 . 0  »  '  .-culs 

•;  t  j  t  e  •  j  3  t  s 

DO  SO  1  =  1. Nr  FT/2 
J  -  I +NFFT  •'  2 
T  £  fi  F  -  C  Y  '  I  > 

C  V  i  I  ■  -  C  Y  <  J  ;• 

1 ; Y < w ) - f  EfiP 
YiJ;  1-1  >*r?CTP 
( •'  I  i  —  {  N  F  F  T  /  2  r  1  ’■  '  f  T  F  0  T  P 

CON r INUE 


.  r  ;•  ■  3  .  y 


■let  *  u  n  2  3  t  e  d  * 


* . 

3'  7 


4 


sir  a  Cat. 
T  V P  E  4  0  C 
c  n  S' m  *,  t  .  • 


Enter  filename  to  save  capstra  d; 


CALL  GETDFN(FIL.NAM) 

Cr’NFIL  (  3  >  F ILNArt  ?  'W  *  >  .ME.  .TRUE.)  GOTO  40 
CALL  PUT DA T • 3  > NFFT  ? TFCTR ?  CY  ) 

or.  sj.  nal  convolved  X  and  Y  data  back  into  3  r  r 


1:3  100  I  - 1  ?  :FFT 

cy  ( i :  =  a  u  x  <  i ; 

Y  ■:  1  *•  -•(  I-l  )  'X TFCTR 
CONTINUE 

VERGE  PROCESSING 
f'-'-E  121 

•  uriflAT  <  /  f  '  Do  you  want  -a  symmetrical  filter  "?  5 
IANS  -  ASK ( ' N  '  ) 

TYPE  133 

T  G  R  M  A  T  (  /  t  '  E  n  t  s  r  bate  l  y  p  a  (1  for*  i  m  ruiss  trsi  n 


F G R M A T  (  '  1  da  t e s  -a y s t e m  i m p u  1  s e  »  2  da t e s  r 3 f e  r e 
A  C  C  E  P  T  #  r  L  0 

•'■IFF  -  <  M  F  F  T  /  2  )  -  1  i  -set  1  imi  t  0  n  t  i  in  e  w  1 

TYPE  1 2 6  > N F F 

F 0 R fl A T  <  /  y  '  Enter  p o -a  1 1 i  v e  t i m e  w  in d o w  (  <  o  r  = 

T  Y  ?  E  :fc  ?  '  A  n  swer  in  ust  b  e  o  n  e  h  3  I  f  o  f  t  o  t  a  1  t  i  in  e  w 
T’-'PE  122 

F  C  R  M  A  T  •'  • j  e  a  1.  red  "or  3  3  y  m  m  eti-ical  filter  t  '  ?  5 

ACCEPT  *f IPCS 

IF \ MFF  . SE.  IPOS)  GOTO  517 

■'■•'PE  *?'  ' 

r  Y  PE  A  '  %  i  J  l  ii  a  elect  e  d  ancae  d  s  p  3  si  ti  v  e  t  i  m  a  p 
TYPE  T  *  '  ceps  t  rum  1  1  : 
fYc,E  •< »  '  ' 


I N  F  3  =  I P  0  S 

IT' IPOS  .EG  . NFF )  INEG  =  IPOS  +  1 
'  -  'IANS  -  N E  .  .TRUE.)  GOTO  5 2 0  !  C h o s e  s y m in et  r  i 

TWPE  131 

l:  0  R  M  A  T  •:  /  7  En  t  e  r  n  amative  t  i  m  a  window  v  I N 
ACCEPT  INEG 
I N  E  G  = I A  8  3 ( INE6? 

IF; INEG  .LE.  NFFT / 2 )  GOTO  520 
TYPE  132 

FORMAT ( '  Width  exceeds  negative  time  point 

in  t  Trt  1  7  >"• 

‘-J  ‘-W  *  'J  \, 


\li 


'MVt  i 


18  !  •-»  r‘  T  =  1 
I  END=  I  F'OS 

DO  160  I=ISTART ? I  £  N  D 

GA TENG=GATENG+CY ( I >*CY <  I  ) 

or < I ) =0 . 

CONTINUE 

I  GTART  =  NF FT-INEG-rl 
I  E  N  I-  -  N  F  F  T 

00  130  1  =  I S  T  A  F;  T  ?  I E M D 

3ATENG-GATENG+CY ( I >*CY( I ) 

C  Y < I ) =0 ♦ 

CONTINUE 

GAI'ENG  =  TOTENG  -  GATENG 
TYPE  1  90 t GW t GATENG 

F ORMA T  (  /  y  '  Enersfy  removed  by  satewidth  of  ’  *  I  3  •  ‘  r- o  : n 
12. 5r‘  joules') 

R  o  t  a  t  e  ■  J  a  t 3 

DO  205  I - 1 ?  N F F T / 2 
J=I+NFFT / 2 
TEMP=CYt I > 

CY( I)=CY< J) 

CY < J)=TEMP 

Y<  J>  =  <  i-i  )*tfctf: 

Y ( I)  =  ( I - <  NFFT/2+1 ) )*TFCTR 
CONTINUE 
TYPE  201 

F 0 R M A T ( / 7  Do  you  wish  to  plot  d 3 1 e d  c e p s  t  r 3  '  y $ ) 

IF  < A 3 K ' ' N ' )  .EG.  .TRUE.)  GOTO  421 
CALL  CLEAR 

CALL  FLUTY CCY^NFFT 70 >  !  RT-11  PLOT 

CALL  P L 0 T X Y  <  Y  *  C Y >  N F F T  »  3 )  i  RT-11  PLOT 

TYPE  420 

F 0 K M A T  {  /  y  '  D o  y ou  wish  to  so v e  -iat e d  cepstr  a  data  '  »  5 
IF  ( A 3 K ( ' N '  )  .EG.  .TRUE.)  GOTO  431 
TYPE  425 

FORMAT ( / r '  Enter  filename  to  save  dated  c  e  p  s  t  r  a  dat 
CALL  GETDFN(FILNAM) 

IF  (0PNFIL(3»FILNAM7 'W' )  . NE.  .TRUE.)  G0T0  422 
CALL  PUTDAT(3»NFFTf TFCTRf CY) 

R  o  t  a  t  e  d  a  t  a 

DO  220  I  -  1  7  N F F T / 2 
J=I+NFFT/2 
TEMP-CY ( I > 

C  Y  (  I )  =  C  Y  (  J  > 

C  Y  <  J  )  -  T  E  M  P 
CONTINUE 
r YPE  *y'  ' 

TYPE  K  7  ...  Now  somputind  I  FT.  Prepare  to  plot  1FT. 


IFUSNX  .  NE,  -1)  GOTO  660 
DU  670  I  =  l.NFFT 
CY< I )  =  -CY ( I  ) 

CONTINUE 

'ert'or  rri  i  nverse  exponential  w  a  i  a  h  t  i  n  2 

IF (LO  . EQ.  2 )  IST=1 

IF ( L  0  .NE.  2)  I3T=IAB3( ISFX ) 

IF  (EW  .EQ.  1.0)  GOTO  526 
F  X  =  1.  . 

DO  235  I  =  I S T ?  N  F  F  T 
CY( I ) -CY ( I ) XFX 
F  X  =  F  X  /  E  W 
CONTINUE 
CONTINUE 

I  3  H  f  T  ~  0 

E  F ( L  0  .EQ.  1)  I SHFT  =  ISFX 
DO  240  I = 1 . NFFT 

Y ( I )=< I-1+ISHFT>#TFCTR 
CONTINUE 

P let  inverse  a 3 1 e d  cerstr u m 


CALL  CLEAR 

CALL  PLOTX.YE  Y.CY.NFFT.3)  !  RT-ll  PLOT 
T Y  p  E  702 

FORMATE/.'  Do  you  want  to  plot  date  on  the 
IF  EASKE'N')  .EQ,  ♦TRUE.)  GOTO  702 
CALL  PLQTXY(Y»CY»NFFTj 1 ) 

TYPE  450 

FORMAT \ / t  '  Do  you  wish  to  s 3 v e  deco n v olve d  data 
IF  E  A  3  K  E  '  N '  )  .EQ.  .TRUE.)  GOTO  461 


FORMATE/.'  Enter  deconvolved  data  filename  t 
CALL  OETDFN E FILNAM ) 

IF  E0PNFILE3.FILNAM. 'W )  .NE.  .TRUE.)  GOTO  45 
CALL  P  U  T  Ei  A  T  E  3  *  N  F  F  T  .  T  F  C  T  R  .  C  Y  ) 

DO  250  1  =  1 ?  N F F T 
CY (  I  )  =  AUXEl) 

CONTINUE 
TYPE  261 

)  U  V  M  A  r  ( .  '  ?  lot  re  m  airdmi  compone  n  t  at  this  a  a 

IF  EASK('Y')  .NE.  .TRUE.)  GOTO  290 

IF  (LO  .EQ.  2)  GOTO  230 

LO  =  2 

GOTO  140 

LO  =  1 

GOTO  140 


•SWTO! 


-  :j r\ iim  i  •»  /  »  :-i n axyze  cepsirs  osis  gi  on  3  01 :  r 

I  F  (ASK  •;  '  Y  '  }  .  NE  *  .  TRUE  .  >  GOTO  300 

T,'PE  330 

FORMAT  <  >  '  ...  Prepare  to  r-  lot  t  o  1 3  i  c  e  -0  s  t 

G  0  T  C  3  4 
TYPE  310 

F  0  R  MAT  (  /  ?  '  !J  i  s  h  to  analyze  another  con  v  o  1  v  e 
I F  : ASK ( ' N ' >  ♦ NE  *  . TRUE . >  GOTO  7 
CALL  EXIT 
£  N  D 


WWW* 


S IJ  B  R  C  U  TINE  F  A  S  T  (  B  *  N  ) 

T h e  subroutines  used  in  this  p  r  o  d  r  a  in  were  1 3 k e n  f  r o m  I 

'  R  roars  ms  For  D  i  d  i  t  a  1  3  i  d  n  a  1  P  rocessind' 

IEEE  Press*  197? 

34  5  East  47  S t  r e 1 1 »  New  York  *  NY  10017 
2  e  o  n sored  by  the  IEEE  Acoustics*  Speech*  and 

Sidnal  Processind  Society 
Lib.  oF  Congress  Cat.  Card  *  79-39023 
IEEE  Book  f  0-37942-128-2  (paperback  ver. ) 

I  0-87942-127-4  (hardback) 

Also  published  by  John  Wiley  S  Sons?  Inc. 

Wiley  0 r d e : ■  3  0-471-05961-7  < paperbac k  v e r . ) 

4  0-471-05962-5  (hardbsc k ) 


DIMENSION  B ( 2 ) 

COMMON  /CONS/PI I ?  P7  *  P7TW0  *  C22  *  S22  *  P  I  2 
PI  I - 4 , *  A  T  A  N  '  1  .  ) 

P  1 8  =  P  1 1  /  8  . 

P  ?  -  1  .  /  S  Q  R  T  (  2  .  ) 

P 7 TWO =2 . #P7 
C  2  2  -  C  0  S  (  P 1 8  ) 

322 =21 N ( PI  3 ) 

PI2=2.*PII 

00  10  1=1*15 
1  - 1 

:l  r  =  2  KY  I 

IF  t N.E0. NT) GO  TO  20 

CONTINUE 

WRITE  (4,9999) 

FORMAT ('  N  --  NOT  A  POWER  OF  2  ') 

STOP 

N4PQW=M/2 

IF ■  M-N4P0W  K2) 40 >  40  *  20 
N  N  =  2 

INT=N/MN 

CALL  FR2TR( INT*B( 1 ) *B( INT  +  1  )  ) 

QU  TO  50 
NN  =  1 


35 


j-'t 


*  L  ■ 


'nvyoViV 


V 


DG  60  IT-1 » N4PQW 
rU'l-H  N*  > 

:  ;  =  n  /  >;  n 

i  •  .1 L  L  F  R  4 1  i  1 T  t  N  N  *  3 '%  1  )  *  S  \  1 N  T  t  1  )  ?  S  (  2  't>  1 N  T  -f  1 

: n r  + 1  •  •  s :  i  ;■ »b<  int+i >  »b<2*int+i  > » b < 3*int+  i  >  > 

CONTI  IJ  £ 


CALL  F  u  R  D 1 ( M ? B  > 
CALL  F  0  R  D  2  (  M  »  B  > 
”  =  B  i  2  ) 

3 (2)- 0 . 

8  <  N  +  1  )  =  T 
8 • N  +  2 ) =0 . 

I1 0  30  1  T  =  4  t  N  t  2 

B ( IT)--B( IT) 

CONTINUE 

RETURN 

END 


SUBROUTINE  F  0  R  D 1  \  M  ?  B  ) 
DIMENSION  B ( 2 ) 

K  =  4 

KL-2 

N=2*#M 

DO  40  J  =  4  ?  N  ?  2 
IF(K-J)20»20»  10 
T  =  B  (  J  ) 

B  (  J  )  =  B  (  K  > 

B(K)  -T 
N-K-2 

IF<K-KL>30»30»40 
K  =  2*J 
KL  -  J 

CONTINUE 

RETURN 

END 


SUBROUTINE  F0RD2(M»B) 

DIMENSION  L(15) »B(2) 

EQUIVALENCE  ( L 1 3  t L ( 1 ) ) j  (Li4»L(2) )  * ( L 1 3  ? L \ 3  )  )  j  (L12?L(4) 
■K  r  (  L 1 0  »  L  (  6  )  )  t  (  L  9  7  L  (  7  )  )  *  (L3»L(3)  )  ?  (  L  7  »  L  (  ?  )  )  *  (  L  6  *  L  (  10)  )  7  (  L  5 
'#(L4»L(12))»(L3fL(13))»(L2»L(14))»<L1»L(15)) 


N=2**M 
L  (  1 )  -N 
DO  10  K  =  2  7  M 
L ( K >=L ( K-l ) / 2 
CONTINUE 
DO  20  K=M,14 
L  (  K  + 1  )  -  2 
CONTINUE 
I  J  =  2 


a 


-  .i 


5  .  *  #  "  v*  -w'  '  «* 


:  2  7  L 1  >  2 


=  J 1  •  L  2  t  L  1 
=J2»L3»L2 
■-J3  t  L4  r  13 
=J4.L5»L4 
:J5»L6»L3 
:  J  6  >  L  7  r  L  6 

*  -J  f  *  L  O  f  L  / 

:JS»L9tL8 
>=J9»L10»L9 
.=J10?L11 ? L10 
>=J11 t L12»L11 
>=J12>L13>L12 
*=J13,L14, L13 


IF  < IJ-JI >30, 120 > 120 
T^B<  IJ-1  ) 

9 •:  u-i  >=B(  Ji-i) 

B( JI-1 )=T 
r=B( i j ) 

B  (  I J  )  -  B  (  J I ) 

B  <  2 1  )  =  T 


r  j  -  j  j  i  2 


RETURN 

END 


SUBROUTINE  FR2TR < INT » BO » 81 > 
DIMENSION  B 0 ( 2 ) *  B 1 ( 2 ) 

DO  10  K= 1 » I NT 
T=BO ( K ) +B1 ( K ) 

B1<K)=B0<  K)-Bl(K) 

B  0  <  K  )  =  T 
CONTINUE 
RETURN 
END 


SUBROUTINE  FR4TR( INT»NN»B0»B1 »B2> B3» B4»B5f B6 » B7 > 

DIMENSION  L<15>»B0(2)»Bl<2)»B2(2)»B3(2)»B4<2)f 
195(2) ?  B  6 ( 2 ) ?  B  7 ( 2  ) 

COMMON  /C0NS/PII»P7»P7TW0»C22»S22»PI2 

EQUIVALENCE  (L15»L(l))f(L14fL(2))»(L13»L(3))»<L12»L<4))»<Lll 
* t < LI  0  J  L ( 6 ) ) f ( L9  f L ( 7 ) ) » ( L8 • L  < 3 ) ) > ( L7  >  L ( 9 ) ) t <  Lo  >  L ( 10 ) ) >  < L5  j  L ( 1 1 > 
*<L4»L<12>)»<L3rL<13)>»<L2»L<14>>f<Ll»L<15>> 


L  '  1  )  -  N  N  /  4 

DC  4  0  K  =  2  r  1  5 

IF  (L  CK-1 > -2 > 10f 20» 30 

L  (  N  - 1  )  -  2 

L  (  K  )  =  2 

GOTO  40 

L ’ K ) - L \ K - 1 )/2 

CONTINUE 


’ 


IV.V.V  V  V  V,‘»*  >'< 


2 ! : 

1  2  0 

2  2-  J  2  *  2 1 »  L  2 

j  0 

1  2 

24-2i»L4>L2 

DU 

1  20 

_  cj  -  J  4  •  L  5  >  .4 

20 

120 

Jc-JjrLfuLD 

DO 

120 

J?-J6*L7  »  L  6 

DO 

120 

J  3  =  J  ^  * L3 ,L7 

no 

12  0 

J  ^  =  J3  >  L9  r  L3 

DO 

120 

J  1  0  -  J  9  «  L  1  0  « 

nn 

12  0 

J  1 1  =  J  1  0  f  L  1  1 

DO 

1  2  0 

J12=Jil»L12 

no 

120 

J13=J12»L13 

n  o 

120 

J14=J13*L14 

O')  120  JTHET--J14»L15>L14 
TH2  -  JTHET-2 
IF(7K2>59»50»90 
tO  60  K= 1 t INI 
TO  =  BO ( K >  +  B2 « K ) 
r l =B1 (K > +  B3 ( K ) 

B  2  (  K  ;  -  B  0  (  K  >  -  B  2  (  K  > 

B.j  (  K  )  =  81  <K  )  -  B3(  K  ) 

£•  0  (  K  }  =  T  0  +  T 1 
Cl  v K ) =  T 0 - T 1 
■2  0 N  r  I N  U  E 

I  F  (  N  N  -  4  )  1 2  0  »  1  2  0  »  7  0 
K0=INT*4+1 
K  L  -  K  0  + 1 N  T  - 1 
DO  80  K  =  K 0  f  K L 
PR=P7*<B1 ( K ) - B 3 ( K  )  ) 
PI=P7*(B1<K)+B3<K> ) 

3  3  (  K  )  =  B  2  (  K  )  +  P I 
31  <K)=PI-B2(K) 

8  2  K  >  =  B  0  (  K  )  -  P  R 
B  0  (  K  )  =  B  0  (  K  )  +  P  R 
CONTINUE 
GOTO  120 

ARG=TH2*PI0VN 

C1=C0S(ARG> 

31=SIN( ARG) 

C2=C1**2-31**2 

32-01*31+01*31 

C3=C1#C2-S1*S2 

S3=C2*S1+C1*S2 

INT4=INT*4 
JQ=JR*INT4+1 
K0=JI*INT4+1 
JLA3T=J0+INT-1 
00  10 0  J=JOfJLAST 
K  - 1\  C)  t  J  J  0 


ngrg 


R  5-31  < J)*S1tB5<K)*C1 
72-=B2(  J>  *C2-B6(SO*S2 
To=B2< J}*S2+B6(K)*C2 
f  3  = B 3 ( j > "f- C 3-8? ( K >  *S3 
T  7  ■=  jf  J  •;  J  )  3  3  r  3  ?  (  K  )  %  C  3 
T  0  =  8  0  f  j  )  t  T  2 
:  !  =  b  •)  t  K  t  1  o 
T  2  -  B  0  \  J  '  --  T  2 
r  -•  =  B  4  :  K  i  -  r  :::• 

:  i  =  R  l  -*• T  3 
7  2  - 1'1  bv  l  . 

1 3  -  R 1  -  T  2 
T  7  =  R  £  -  T  7 
30 ( J)=T0+T1 
0  7 (K)  =  M  +  T5 
■  6 <  K ) -T0-T 1 
B1 ( J)=T5-T4 
32<J>=T2-T7 
8 3  < K >=T6+T3 
P4  < K ) =72  +  T7 
3  3  (  J  >  =  T  3  -  T  6 
CONTINUE 
J  K J  R  y  2 
J I  =  J I  -2 

IF< JI-JL) 110t 110» 120 

JI=2*JR-1 

JL  =  JR 

CONTINUE 

RETURN 

END 


SUBROUTINE  F3ST<B?N) 

DIMENSION  B ( 2 ) 

COMMON  /CONST/P  1 1 »  P7  >  P7TU0  » C22 
PII=4.#ATAN< 1 . ) 

P I  8  =  P 1 1  /  8  . 

P  7  =  1 ./SORT ( 2 . ) 

P7TW0=2.*P7 
C  2  2  =  C  0  S  <  P  I  3  ) 

3  2  2  =  8 1 N  (  P 1 8  ) 


12=2. *P I 


DO  10  1  =  1 » 15 
M  =  I 

N  T -2%% I 

IF (N.EQ.NT)OO  TO  20 

CONTINUE 

!J RITE  <  4,999?) 

FORMAT ( '  N  --  NOT  A  POWER  OF  2 
STOP 

B  •:  2  )  =  6  (  N  + 1 ) 

DO  30  I  =  4  »  N  j  2 
B  < I )  =  -B ( I  ) 

CONTINUE 

DO  40  1=1? N 

B  <  I  >  =  B  (  I ) /FLOAT  '  N  ) 


m 

i 


i#« 

tJ'M 


CONTINUE 


N 4PDW=M/2 


CALL  F 0 R D 2  Li  j  6  ) 

C  A  L  L  F  0  R  0 1  (  M  ?  3  ) 

IF  (  N4P0W  .  EQ  .  0  )  GOTO  60 
NN -4 NN 

DG  50  IT=1»N4P0W 
N  N  ~  N  N  /  4 
INI =  N  /  N  N 

CALL  FR43YN(INTfNN>B<l)»B(INT+l)»B<2*INT+l} 
1  t  B<  3*INTrl ) j  3 « 1 ) >  S( INT  +  1 ) >BC  2*INT-!-l ) 

1  ? 3 < 3*INT-hl )  ) 

CONTINUE 


IF  <  M  -  N  4  P  0  W  *  2  )  8  0  •  3  0  »  70 
INT=N/2 

CALL  FR2TR< INT»B( 1 ) »B< INT+1) ) 

RETURN 

END 


SUBROUTINE  FR4SYN( INT »NN»B0»B1 >B2»33»B4 »B5 >Bo » B7 > 

DIMENSION  L< 15)  »B0(2) fBI (2) » B2(2) rB3(2) r B4<2>  7  35 (2) 7  36 ( 2>  r 37(2) 
COMMON  /CONST/PI  I » F7 » P7TW0 » C22 » S22 » PI2 

EQUIVALENCE  (L15fL<l))f(L14»L(2))»<L13»L<3>)f<L12»L<4))»<LllfL(5) 
#.  (L10»L(6) ) » (L?»L(7) ) » (L8»L<8) ) » (L7»L<9> ) » (L6»L( 10? ) »  <L5»L( 11 ) ) » 
*(L4»L<12>)f{L3»L(13))f(L2»L(14)),(Ll»L(l5)) 


L  (  1  )  -  N  N  /  4 


DO  40 


IF  (L(K-l)-2)10t20»30 
L ( K-l ) =2 


GOTO  40 

L  fO  =  L  (  K  -  1  >/2 

CONTINUE 


R 1 0  V  N  -  F’  1 1  /  F  L  0  A  T  (  N  N  ) 
JI--3 


JR  =  2 


2  7  L 1  7  2 

J 1 >  L2  >  L 1 
J2  ?  L3  7  L2 
J3  »  L4 , L 3 
J  4  »  L  5  f  L  4 
J5  7  L6  7  L5 
J6  7  L 7  *  L6 
J7»L3rL7 
J8»L?»L8 
)=J?»L10f L? 

i  =  J 1 0  »  L 1 1  >  L 1 0 

:-JllfL12»Lll 

>=J12»L137L12 


